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Summary 

An approximate axisyrnmetric method lias been developed which can reliably 
calculate fully viscous hypersonic flows over blunt-nosed bodies. By substituting 
Maslen’s second order pressure expression for the normal momentum equation, a sim- 
plified form of the viscous shock layer (VSL) equations is obtained. This approach can 
solve both the subsonic and supersonic regions of the shock layer without a starting so- 
lution for the shock shape. The approach is applicable to perfect gas, equilibrium, and 
nonequilibrium flowfields. Since the method is fully viscous, the problems associated 
with coupling a boundary-layer solution with an inviscid-layer solution are avoided. 
This procedure is significantly faster than the parabolized Navier-Stokes (PNS) or 
VSL solvers and would be useful in a preliminary design environment. Problems 
associated with a previously developed approximate VSL technique are addressed 
before extending the method to nonequilibrium calculations. Perfect gas (laminar 
and turbulent), equilibrium, and nonequilibriuin solutions have been generated for 
air flows over several analytic body shapes. Surface heat transfer, skin friction, and 
pressure predictions are comparable to VSL results. In addition, computed heating 
rates are in good agreement with experimental data. The present technique generates 
its own shock shape as part of its solution, and therefore could be used to provide 
more accurate initial shock shapes for higher-order procedures which require starting 
solutions. 
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1 Introduction 

Ongoing investigations into configurations such as the Aeroassist Space Transfer 
Vehicle (ASTV), the National Aero-Space Plane (NASP), and the Personnel Launch 
System (PLS), with their associated high altitude, high speed environments, have 
sparked renewed interest in hypersonic aerodynamics [1]. Since extensive computer 
run times prevent more exact approaches from being used in the preliminary de- 
sign environment, there is a continued interest in developing improved engineering 
methods. 

For large Reynolds numbers, the shock layer consists of a large inviscid region 
and a thin boundary layer near the body. Typically, an inviscid solution to this 
outer region is coupled with a boundary-layer technique (e.g., Ref. [2]). For flows 
about blunt-nose bodies, the shock slope changes rapidly in the streamwise direction 
creating strong entropy gradients in the inviscid flow which must be accounted for 
when coupling the inviscid-layer and boundary-layer solutions. 

The viscous region encompasses a significant portion of the shock layer at the lower 
Reynolds numbers encountered by hypersonic vehicles at high altitude. Most inviscid- 
layer /boundary-layer approaches neglect the boundary- layer effect on the outer layer, 
so in some cases a simple coupling technique may not give satisfactory results. Ac- 
curate predictions of the flowfield properties may be obtained from solutions to the 
Navier-Stokes [3], the parabolized Navier-Stokes [4, 5, 6], or the viscous shock layer 
[7, 8, 9, 10] equations. The full Navier-Stokes equations are typically solved us- 
ing a time-marching procedure in order to properly model their elliptic behavior (a 
very costly computation). The parabolized Navier-Stokes (PNS) and viscous shock 
layer (VSL) equations are derived from the steady compressible Navier-Stokes equa- 
tions. Both equation sets have been parabolized in the streamwise direction so that 
a solution may be advanced downstream using spatial marching techniques. These 
formulations account, among other things, for the nonzero normal pressure gradients 
which are neglected in a classical boundary-layer approach. Unfortunately, the com- 
putational requirements of existing methods for solving these equations exceed that 
which can be tolerated in preliminary parametric design studies. 

For high energy flows, the assumption of perfect gas behavior is not valid since 
chemical reactions occur. Using the perfect gas assumption results in much higher 
post-shock temperature predictions than are realized in the flight environment. This 
reduced temperature is due to energy being absorbed by endothermic chemical reac- 
tions. The classical approach is to assume that these reactions occur at a fast enough 
rate, relative to the time scales of the flow, that an equilibrium composition has been 
achieved everywhere in the shock layer [8, 1 1]. 

Chemical reactions are a result of sufficiently high energy molecular collisions 
which, in air for example, cause the dissociation of O 2 and N 2 into 0 and N , along 
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with the production of NO. Such reactions require a finite amount of time to occur. 
Thus, a finite amount of time is required for the related chemical changes to take place. 
For equilibrium flows, the assumption is made that the residence time of the fluid in 
the shock layer is large relative to the time scales of the chemical reactions. Hence, 
the reactions occur nearly instantaneously. The assumption of equilibrium flow is 
not correct for many cases of interest, particularly with reentry trajectories. This is 
evident from the Catalytic Surface Experiment [12] flight data, which illustrate the 
nonequilibrium phenomenon via the rather noncatalytic nature of the Space Shuttle 
surface. Thus, consideration of these nonequilibrium effects in the design process 
is desirable. VSL solutions with finite-rate chemistry have been documented in a 
number of references [8, 13, 14, 15, 16, 17]. Unfortunately, the restrictive central 
processing unit (CPU) requirements cited earlier are intensified when nonequilibrium 
effects are included. Therefore, a more approximate approach which accounts for the 
effects of finite-rate chemistry would be useful. 

In 1964, Maslen [18] published a simple inverse method (shock shape prescribed, 
body shape to be determined) for calculating the inviscid flowfield within the shock 
layer surrounding a smooth axisymmetric body. The direct problem, which is the 
more straightforward application, requires iterating the shock wave shape until the 
desired body shape is obtained. A closed-form expression for the local pressure in 
terms of the stream function is obtained by approximating the integral of the normal 
momentum equation. With the original approximations, the effects of the velocity 
normal to the shock are neglected. Using this first-order expression, Grose [19, 20] 
developed an inviscid, nonequilibrium method for calculating flows in Earth, Martian, 
and Venusian atmospheres. 

In a later effort, Maslen approximated the normal velocity’s contribution to the 
pressure relation to obtain a second order approximate integral of the normal mo- 
mentum equation [21]. This expression better accounts for the recompression across 
the shock layer. Zoby and Graves [22] coupled this improved relation with the itera- 
tive scheme of Jackson [23] to solve the direct problem for a variety of blunt bodies. 
Later this approach was extended [21, 25] to include an approximate technique for 
calculating heating rates for bodies at angle of attack. Ref. [26] provides a review of 
additional approximate heating methods. 

In a more recent work, Grantz and DeJarnette [27, 28, 29] employ the second 
order Maslen pressure relation, along with a simple linear expression for the normal 
component of velocity, in an approximate VSL approach. Boundary-layer-like viscous 
terms are added to the inviscid streamwise momentum and energy equations to obtain 
a parabolic equation set analogous to the full VSL equations. The viscous terms 
are retained across the entire shock layer, thus avoiding the problems which can be 
encountered in the coupled inviscid-layer/boundary-layer approach. Furthermore, 
since the shock shape is part of the solution, no initial shock shape or smoothing of 
intermediate shock shapes (both required for the VSL technique) is necessary. 

The method of Ref. [29] did not fully yield the anticipated reductions in CPU 


3 


requirements as compared with the full VSL equations. Based on a review of the 
algorithm, several areas of concern in the original technique were recognized: 1) too 
many iterations are required for shock shape convergence; 2) the solution is incon- 
sistent near the stagnation line; 3) the grid-point spacing across the shock layer may 
yield oscillations in the shock layer property profiles. Thus, it became apparent that 
a different approximate VSL approach should be explored. 

The first portion of this paper discusses that effort, and details the differences 
between Ref. [29] and the new approach [30]. Three major differences between these 
techniques are: 1) more efficient algorithms than those employed in Ref. [29] are used 
to generate the shock shape; 2) the governing equations give a consistent limiting 
form on the stagnation line; 3) and the spacing across the layer is related to physical 
distance rather than the stream function, as in Ref. [29]. Comparisons of results and 
run times are made using the current approach, the method of Ref. [29], and a full 
VSL solver [9]. Comparisons with experiment are also made. 

The second part of this document concerns nonequilibrium flowfield calculations. 
A seven-species finite-rate chemistry model for air [15] has been incorporated in the 
approximate VSL solver discussed above. Comparisons of results and run times axe 
made using the current approach and a full VSL solver [14], Comparisons with Space 
Shuttle flight data [31] are also made. 
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2 General Analysis 

Ref. [29] discusses the development of an approximate viscous shock layer tech- 
nique which is solved in a shock-normal coordinate system. The approach begins with 
an approximate inviscid flowfield solution based on the second-order pressure relation 
of Maslen [18], including an approximate representation for the normal component of 
velocity. In order to approximate a full VSL solution, boundary-layer-like terms are 
added to the inviscid energy and streamwise momentum equations. 

Since the goal of the present work is to develop an approximate VSL algorithm, 
the current approach starts with the full VSL equations written in orthogonal curvi- 
linear coordinates. First, the equations are cast in a shock-normal (rather than the 
traditional body-normal) coordinate system. This step is necessary in order to fa- 
cilitate the use of Maslen’s pressure relation which replaces the normal momentum 
equation. Since the remaining equations are unchanged, the normal component of 
velocity is found by solving the continuity equation rather than assuming a profile 
as was done in Ref. [29]. If the form of the energy and streamwise momentum equa- 
tions used in Ref. [29] are compared to the full VSL equations in the shock-normal 
system, it can be noted that several “higher-order” terms (involving the variation of 
the metrics across the shock layer) have been neglected. These higher-order terms 
are retained in the present study. 

The governing equations (along with their stagnation line formulations), bound- 
ary conditions, and surface quantities are presented in this chapter. In addition, key 
differences between the current approach and Ref. [29] are cited. The subsequent 
four chapters deal with applying these equations to perfect gas, equilibrium, turbu- 
lent, and nonequilibrium flows, respectively. Unless otherwise noted, the governing 
equations as presented in this chapter are employed in those applications. As a final 
comment, dimensionless variables (see Appendix A) are employed throughout this 
paper. Dimensional quantities are denoted by a superscript *, except in the results 
chapter where this superscript is omitted. 


2.1 Viscous Shock Layer Equations 

The Navier-Stokes equations written in curvilinear coordinates for axisymmetric 
flow are presented in Appendix B. If v and n are assumed to be order e (where t 
is the Reynolds number parameter), and terms greater than order e are neglected, 
the VSL equations [7] are obtained (see Appendix B for details). They are presented 
below, excluding the normal momentum equation, in nondimensional form. 
continuity: 

jT (puh 3 ) 4- {pvhihz) = 0 


( 2 . 1 . 1 ) 
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s - momentum 


' u du du uv dhi 

_J_ '» i — -l- 

Hi os on n i on 


energy: 



( 2 . 1 . 2 ) 


u dhi 
h\ dn 


hi dn 
du u 

1 dn dn ) \hi dn h, 
where the Reynolds number parameter is defined as 


(2.1.3) 


2 _ 

Poo^oo^noae 

and Pr is the Prandtl number. The dimensional reference conditions are given in 
Appendix A. 

At this stage, this equation set could be transformed to either a body-normal or 
shock- normal system. The shock-oriented coordinate system is shown in Figure 2.1. 
Metrics associated with the shock-oriented coordinate system (where s is the distance 
along the shock ) used in this effort are 


hi = 1 — nn s 


h 3 = r = r s — n cos 1% (2.1.4) 


where n is the inward normal distance 
By definition, 


from the shock and 

_d£ ± 

ds 


k s is the shock curvature. 


(2.1.5) 


Note: In addition to Eqs. (2.1.1) through (2.1.3), the standard VSL 

equations employ the normal momentum equation: 


/ u dv do u 2 dhi \ dp 

^ \hi ds 1 dn hi dn ) ^ dn 


( 2 . 1 . 6 ) 


This set of equations is solved in a body-oriented coordinate system (where 
s is the distance along the body). The metrics are defined as 


hi = 1 + nn b h 3 = r = n + n cos Tt 

where n is the outward normal distance from the body and Kb is the body 
curvature (which, by definition, is —dTbjds). 
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Figure 2.1. Flowfield geometry and curvilinear coordinate system. 


In order to facilitate the solution of the governing equations, a transformation 
to normalized coordinates is performed to yield a constant number of points across 
the layer from station to station. Using the shock properties, define the following 
computational variables: 


f = * 


u = 



n 

nb 


where n b and the shock properties are functions of ( only. From its definition, r, n = 0 
on the body and r, n = 1 at the shock. The chain rule of differentiation gives 


Since 


d dZ d dr] n d 

ds ds dZ ds dr) n 

d 

dn 

I- 

dr, n 

ds 

I- 

dr,» 

On 


a * iJL 

dnd£ dn dr, n 


n drib 

nl ds 
1 


n b 


the derivatives in this coordinate system are 


d_ 

ds 


d 

dt 


j] n - 1 dii b 


d 


dn 


i d 

rib dr] n 


n b dZ drj n 


(2.1.7) 
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In the transformed coordinate system employed here, the metrics, Ecj. (2.1.4), become 
h x = 1 + n 6 (r;„ - 1)«. h 3 = r 3 + n b {ri n -1) cost, (2.1.8) 

Their partial derivatives in the normal direction are 

(2.1.9) 


ah _ _ _ _2_oh 

dn 3 n b dr] n 


0h 3 „ 1 dh 3 

= — cos 1, = — 


On 


n b drj n 


The governing equations in this normalized shock-oriented coordinate system can now 
be written as 
continuity: 


I wa) - <**•> - *£ <**•*•> - 0 


( 2 . 1 . 10 ) 


£ - momentum: 


P*P 


. h\ 


Oil 

'ae 


du 3 


u 3 — + w-rr - 


rj n — 1 </«6 du n 

rift d£ dr/ n rift 


dti u 5/ii 

+ 


dr) n hi dfjn 


1 

1 

dp 

T] n - 1 drift 3p 

_ e 2 

( 3 f / u 5/ii\ 

/?> J Titj 

dl~ 

rift d£ <9r/„ 


( drjn [ \dp n h x dr) n J 


energy: 


+ P 


PsP 


2 dh x 1 dh 3 

+ 


.h x drj n h 3 di)r 


( U 3 U 

'dh 

Tj n — 1 drift dh 

\h x 

[d( 

rift d£ drj n 


u a u 

h i 


dp T] n — 1 drift dp 

di ' 


du u dh x \ ) 

drjn h\ drjn J J 

v dh\ 
n b dr) n ) 

v dp 


(2.1.11) 


rife d£ 5r/„J rift drj n 


+ 


( 2 . 1 . 12 ) 


5 


P 


dh 


+ 


p dh 


l dhi 1 5*3 \ 2 / du u dhA 2 ) 

hi drjn + h 3 dr] n ) + * \d»?n h x dr] n ) J 


jdr/ n [PrcfynJ Prdr; n \/ii 

Equations Eqs. (2.1.11) and (2.1.12) can be cast in the following standard form 
for parabolic partial differential equations: 


d 2 W dW . dW „ 

4o^T7 + + A t W + A 3 + A 4 — = 0 


(2.1.13) 

o o t - o i i ' V f 

dpi OTj n 

where VE represents the dependent variables u and /i, respectively. The values for the 
coefficients are listed below: 

£ - momentum: 

An = --5/1 
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At = — , 


nt 


v u,u Tj n — 1 drib 

iib h\ Tib d£ 


A 2 — — o 
nr 


' 1 dhi 
h\ dr) n 


t , / 1 1 < 9 /< 3 \ 

dp f 1 dhi 1 dh 3 \ p a p fdu, v dh\ 
prj n + ^ \hi drj n + h 3 drj n ) + hi \ d£ n b dij nj 

«i U C'c J ori n \u 3 J 


( 2 . 1 . 14 ) 


A P U 

A 4 = p a u a — 

h\ 


energy: 


>1, = 


~ ,2 


e 2 p 


n l Pr 


c 2 

' d 

{ P \ 

, P i 

(l dhi 1 dh 3 V 


V ( u a u rj n 1 drib 

n l 

dr) n 

\ Pr) 

+ Pr 1 

[hi dy n h 3 dr) n J 

-PsP 

rib hi rib d£ 


Ad — 0 


( 2 . 1 . 15 ) 


A3 = — 


■u a u 

~ih 


Op 

Tjn - 1 drib dp 

, v dp 2 u] 

fdu 

u dhi \ 


rib d£ dg n 

+ * t P 2 

n b dr) n nl 

\dr/ n 

hi dr) n J 


A\ = p s u , 


pu 

% 


2.2 Maslen’s Pressure Equation 


As mentioned in the previous section, the normal momentum equation is solved as 
part of the standard VSL solution. In the approximate VSL method, Maslen’s second- 
order pressure equation is used in lieu of numerically integrating the n-momentum 
equation. The derivation of this approximate closed-form expression is presented in 
Appendix C. The result is repeated here: 


,, , . Hsr a u a v a sin 

Pit, T ))=P> J r — Tj— (»? - 1) 7 




where 

and is the stream function. 

Maslen’s pressure expression offers an attractive alternative to numerically in- 
tegrating the normal momentum equation to solve for p, since it expresses p as a 
function of the shock properties at the current station and t] only. This relation is 
very simple as compared with the normal momentum equation which it supplants. 
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The variable ?/ should not be confused with the transformation variable rj n defined pre- 
viously. The relationship between those distinct quantities is derived in Appendix D, 
and is repeated below: 


I* rr)n 1 f 

'i hv~ n l P* u l pu{T}n-\)di] n +n b p,u,r, I pudr] n (2.2.2) 

L Jo Jo 


Maslen’s relation fosters analytic relations for the partial derivatives of pressure 
(which appear in the governing equations) as well. Differentiating Eq. (2.2.1) with 
respect to rj gives 

dp _ KsTjU, _ i>s sin r, T Kjr, 

dp ~ 2 2 1 cos T, 

Differentiating Eq. (2.2.1) with respect to £ gives 


(2.2.3) 


dp 

ft 


dp, 

d£ 


+ juj/Cs sin + u,r + k]t , sin Taj 


7/ 2 — 1 

f v, sin T s 

+ k 3 sin T, ( 


4 

[ COS I 5 

d£ \ 

COS I j / 


+ 


(l + -^-) 

V cos 1 , / 


sin T 


dv, 

'ft 


V,K s CO S r, | 


(2.2.4) 


As noted, Eq. (2.2.4) is the partial derivative of p along lines of constant rj. The 
following expression relates this quantity to the partial derivative of p along lines of 
constant p n (see Appendix C): 


where 


and 


Further, 


dp 

ft 


Vn 


dp 

ft 


drj dp p n ~ 1 drib dp 
+ TTTT + 


ds drj rib d^ drjfi 
dr/ p,pu,uh 3 dn I 


ds 


di 


dn 

ft 


Q-v) ^ 

2 cos r s 


1 - 


K„r. 


cos T.J 


dp p,pu a uh 3 rib dp 


dt) n 




drj 


(2.2.5) 

( 2 . 2 . 6 ) 

(2.2.7) 

( 2 . 2 . 8 ) 


2.3 Wall Boundary Conditions 

The no-slip condition is applied at the wall so that* 

u w = v w = 0 

In addition, the wall temperature (T w ), which can vary along the wall, must be 
specified. 
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2.4 Stagnation Line 

The stagnation line is a singularity in the standard VSL equations. Typically, this 
singularity is handled by representing the dependent variables on the stagnation line 
with a truncated power series of £. Substituting these expressions into the governing 
equations yields a set of ordinary differential equations for the stagnation line. How- 
ever, since truncated expansions were used in their formulation, these stagnation line 
equations are not strictly consistent with the general governing equations. 

In the shock-oriented system, an explicit limiting form of the governing equations 
as £ — + 0 can be calculated, since Maslen’s relation provides analytic derivatives 
of the pressure. Thus, a set of ordinary differential equations which is consistent 
with the general governing equations can be obtained. The form of the streamwise 
momentum equation used in Ref. [29] contains a group of terms which vanish on 
the stagnation line provided the normal momentum equation is satisfied. However, 
Maslen’s expression does not satisfy the normal momentum equation, so a correction 
term must be carried along throughout the shock layer in order to force the solution to 
smoothly approach its stagnation value. Despite these efforts, the stagnation solution 
was not consistent with the remainder of the layer. In the current research, the 
streamwise momentum equation does not contain this group of terms, so no correction 
terms are required. Results of Chapter 9 show that with the current approach, the 
problem with the limiting form of the equations as they near the stagnation line is 
successfully addressed. The governing equations for the stagnation line are presented 
below. 


2.4.1 Viscous Shock Layer Equations 

The limiting values of the geometric quantities are given in Appendix E. First 
consider the continuity equation. Without loss of generality, Eq. (2.1.10) can be 
rewritten as 


^ ( p 3 u a n b puh 3 ) - [(Vn - 1 ) piih-i) - p (pvh x h 3 ) = 0 (2.4.1) 

On the stagnation line, 


{p,n b pulh) + ^puh x ^{u s r s ) - PsU^^-Kr ) „ - l)puh x ] 


di 

~ Ps "din = 0 


dt d Vn 


so that 


p r uh x lim < /V4 


u s dr s du s 

r s d( + d( 


(2.4.2) 


} ~ Pa °Wn PVh = 0 (2 ' 4 ' 3) 
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Further applying these limits to Eqs. (2.1.11) and (2.1.12) yields the following set 
of equations which are valid along the stagnation line: 
continuity: 

2>is 0 nb 0 ps 0 puhi - p a 0 ^- (pvhfj = 0 (2.4.4) 


£ - momentum: 

P< oP 


( u 2 du 3 i 

d£ n bQ [dr} n ' hi dy n \ 


Du u dhi 

+ 


\ 1 f 1 \djp rjn - 1 dnb dp 1 1 

/ hi n b d£ dr) n \] 

2 f d T / du_ _ (&1L _ (o 4 k\ 

o 2 L \% hidTj n )\ hidr) n \dr) n hidTj n )J 


energy: 


Ps 0 P v 


dh 


+ V 


dp 


dT) n drjn 


e 

n b n 


/ « 

p dh 

\ df) n 

Prdrf n 


+ 


2 p dh 


hi Prdr ) , 


q d/n 

'n / 


It follows that Eqs. (2.1.14) and (2.1.15) become 
£ - momentum: 


Ai — 


/ 4i = - 

1 d/q 


j 4 q — 

dp 


+ 2 


2 ^ 

n 6o 2 

p d/q 


« 6 n 


As — 


Aj d? 7 „ 

1 


« 6 0 l 
dp 


+ 2 


drjn hi drj 
p dhi 


~ P*oP 


«fc n 


dt] n hi dl) n 


h i 


d£ J n 6(J dr) n €-o ( u 4 d£ J 


Al = 0 


(2.4.6) 


< 2 - 4 - 7 > 


energy: 


Ao — — 


ilii 

«t 0 Fr 


^ = 

n bo 


d ( p \ 2 p dhi 


TT + 


dr) n \PrJ h { Prdj] n \ 
A 2 = 0 
dp 

dp n 

Ai - 0 


^ 0 ^ 


(2.4.8) 
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2.4.2 Maslen’s Pressure Equation 

On the stagnation line, equation Eq. (2.2.1) reduces to 

7) = p> 0 ~y (v 2 ~ l) 

and Eq. (2.2.2) becomes 
1 


IV = n b n 


From Eq. (2.2.3) 

Further, Eq. (2.2.4) becomes 

limU| 

<-*0 U, d( 


p3 ° K, ° 2 Jo PU ~ ^ (iT}n ] + ” 6 0 b -0 J 0 P ud7 lri^ 

Op 


dp 


= ~*v/ 




Utilizing results from Appendix E, this can be written as 

HmU g 

I u a d£ 


7/ 2 — 1 


where the terms 


Bn./!* 

[ u, d£ 




(2.4.9) 


(2.4.10) 


(2.4.11) 


(2.4.12) 


(2.4.13) 


are defined elsewhere. 

Continuing, from Eq. (2.2.5) the limiting form of the streamwise pressure gradient 
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where 



= 2 p, Q puh x lim 



Taking the limit of the streamwise derivative of n gives 

r cos T, - r a 



= ( 1 -,)^i iin j cos V- — } 

2 <-»o [ u*r s J 


so that (see Appendix £1) 


lim 




(2.4.15) 


(2.4.16) 


Further, the derivative of p with respect to r/„ is 
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3 Perfect Gas Analysis 

The governing equations given in Chapter 2 must be supplemented by the equa- 
tion of state. In addition, the thermodynamic and transport properties are required 
for closure of the equations set. This chapter supplies those quantities for perfect gas 
flows. In addition, the perfect gas shock jump conditions are given in Section 3.3. 
Perfect gas results generated using these governing equations are presented in Chap- 
ter 9. 


3.1 State Equation 


In its general form, the equation of state is 

P* = P* 


For perfect gas flows, 

p* = p*il*T* (3.1.1) 

where for standard air ( and 0-i only), TZ = R* air . Written in nondimensional form 
(see Appendix A), 

P = pRaxrT = p-~—T (3.1.2) 


where 7 is the ratio of specific heats. 


3.2 Thermodynamic and Transport Properties 

For perfect gas flows, the specific heat has a constant value: 

NR* 

c; = ~ (3.2.1) 

where in nondimensional form (see Appendix A) this becomes 

<V = yrf = 1 ( 3 - 2 - 2 ) 

As a result, the definition of enthalpy for perfect gas flows is simply 

h* = C;T * t (3.2.3) 

h = C P T = T 


so that nondimensionally 


(3.2.4) 
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Sutherland’s law [32] is employed to calculate the viscosity: 


p * = 2.27 x nr 8 - r 


T* + 198.6 


[slug/ft- s] 


(3.2.5) 


where T* is in degrees R. The Prandtl number ( Pr ) is the ratio of viscous diffusivity 
to thermal diffusivity: 

Pr= ^ (3.2.6) 

a! 

For a perfect gas, this quantity is constant throughout the shock layer, and in this 
investigation is assumed to be Pr = .72. With C p , /i, and Pr as defined above, the 
thermal conductivity, k, can be determined from Eq. (3.2.6). 


3.3 Shock Properties 

The shock jump conditions are given by the Rankine-Hugoniot relations. For 
the special case of perfect gas flow, these expressions take the following form (see 
Appendix F for details): 


2 sin 2 F, 7 — 1 

P ’ 7+1 7(7+1) Ml, 

(3.3.1) 

(7+ 1) Ml sin 2 T. 

Ps 2 + (7 - t)Af£, sin 2 T, 

(3.3.2) 

rp 7 

1 

ps 7 - 1 

(3.3.3) 

sin r a 

V s ~ 

(3.3.4) 


The streamwise derivative of u, appears in the ^-momentum equation. Further, 
the derivatives of p s and v s with respect to £ are required for Maslen’s analytic 
expression for dp 3 /d£. These gradients are given by (see Appendix F) 


du s 

dt. 


dv s r 

— = k 3 cos 1 s 

d£ 


*3 sin F 3 
4 


dp, 

d{ ~ ( 7 +l) 


1 

[(7 + 1 ) sin 2 r , p,\ 
\ 


k s sin r, cos r s 


(3.3.5) 

(3.3.6) 

(3.3.7) 


On the stagnation line, the shock jump conditions take the form of the normal 
shock relations: 

7 - 1 


V ’° 7 + 1 7(7+1) Ml 


(3.3.8) 
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t} 


(l+HjC 

2 + h-i)Ml 


T. 


°o 


PfQ 7 
Ph> 7 ~ 1 


= 


Avo 


In addition, the stagnation line limiting forms of the shock gradients are 



(3.3.9) 

(3.3.10) 

(3.3.11) 

(3.3.12) 

(3.3.13) 

(3.3.14) 


3.4 Surface Quantities 


In the design environment, estimates of the lift, drag, and heating rates are de- 
sired. Thus, a design tool should provide values for body pressure, skin friction, and 
heat transfer rate. In this effort, the body pressure is supplied by Maslen’s relation. 
Relations for the skin friction and heating rates are given below. Gradients along the 
shock-normal lines are used in lieu of body-normal values. 

The skin friction coefficient is defined as 


9 



where the shear stress at the wall is 


r = it 

tv r u 



In terms of the nondimensional variables (see Appendix A), 


(3.4.1) 


(3.4.2) 


so that 




V* 

r oo 


R: 


) 


Fw 

iib 




du 

Oljn 


\w 


(3.4.3) 


(3.4.4) 
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The heat transfer rate (see Appendix C) to the body is 

dT* 

n* = —k* 

w dn* 

W 

so that in nondimensional form, 

2 dT _t\ dT 

€[w € K W - foyj 

dn w n b d Tj n w 

Since 

dh. = C p dT 

this may be written as 

( 1 k w Oh 

Hu> — ' 7T "o 

lib Op CyTJ n 

or 

_ c 2 i-i u , dh 

<lw n b Pr w dr) n 


W 


( 3 . 4 . 5 ) 

( 3 . 4 . 6 ) 

( 3 . 4 . 7 ) 

( 3 . 4 . 8 ) 
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4 Equilibrium Analysis 

The governing equations as presented in Chapter 2 are readily applicable to equi- 
librium flowfields. In fact, the form of the continuity equation and Maslen’s equation 
remains unchanged for all flow regimes considered in this paper. The state equation, 
the thermodynamic and transport properties, and shock jump conditions for equilib- 
rium flows are provided in this chapter. Chapter 9 presents equilibrium results based 
on these governing equations. 


4.1 State Equation 

For a nonreacting gas, R is constant since M is constant. However, when 
chemical reactions become important, M varies within the fluid, and using R = R air 
in Eq. (3.1.2) is no longer valid. There are a variety of approaches for handling 
equilibrium flows [33, 34] which account for the variation of R . In this study, 
Eq. (3.1.2) is modified to account for the effects of chemical reactions: 

P = pZRairT (4.1.1) 

The term Z is the compressibility factor: 

ry Af at > 

This factor is evaluated through a table look-up procedure of values calculated using 
Hansen’s expressions [35]. 


4.2 Thermodynamic and Transport Properties 

For equilibrium flows, both the Prandtl number and specific heat of the fluid vary 
within the shock layer. In addition, Sutherland’s viscosity relation is no longer valid. 
In past efforts (see [36], for example), various techniques [37, 34, 38, 39, 40] have been 
employed to define the thermodynamic and transport properties of equilibrium flows. 
The energy equation as written in Chapter 2 requires that the equilibrium Prandtl 
number be used (see Appendix G). Hansen’s paper [35] includes expressions for h , 7 1 , 
/i, and Pr eq . This is the model used to define thermodynamic and transport properties 
for the solutions presented in this paper. As in the calculation of Z, these relations 
are used to generate a table of data which spans a large range of values for p, ft, 
and T. These tabulated values are then interpolated as needed during the numerical 
solution of the governing equations. At the body surface (where T w is known), the 
table look-up procedure takes the form 
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Hw 


h W [Pw ? -1 w ) 
W ( pw ? ^ If ) 


— Z w (p w , T w ) 

P r eq w ~ P r eq w {PwiT w ) 


Elsewhere, since h is known, the implementation is of the form 


T = T(p, h) Z = Z(p,h) 

P = p(p, h) Pr eq = Pr eq (p, h ) 


4.3 Shock Properties 


With perfect gas Hows, closed-form relations can be used to determine the values 
of the jump conditions. However, no such simple expressions can be formulated for 
chemically reacting flows. Therefore, the following set of equations must be solved in 
an iterative fashion (see Appendix F for details): 


Sill I , 

V a - 

(4.3.1) 

ps 

Ps = Poo + sin 2 r s ( i - — ) 

(4.3.2) 

V pJ 

sin 2 F s ( 1 \ 

h, - /loo + 2 y p2j 

(4.3.3) 

T a {p a ,h a ) Z a = Z a (p s , h a ) 

(4.3.4) 

Ps 

P — 

7 1? .T 

(4.3.5) 


The iterative procedure is as follows: 


1) Begin by setting p s = 10. 

2) Use Eq. (4.3.2) to calculate p s . 

3) Find h a from Eq. (4.3.3). 

4) For these values of p a and h a , find T„ and Z a from Hansen’s model. 

5) Determine a new value of p s from Eq. (4.3.5). 

6) Steps (2) through (5) are repeated until the relative error between 
the new and old values for p 3 is less than the prescribed tolerance. 

7) After convergence, compute v a from Eq. (4.3.1). 


The derivatives of u a , p 3 , and v a with respect to £ are also required. The definition 
of the streamwise derivative of u a as given in Section 3.3 is valid for equilibrium flows. 
F'inding the gradients of p a and v 3 , however, presents a challenge. Fortunately, these 
gradients only appear in Maslen’s streamwise pressure gradient. As will be discussed 
later, a simple two-point backward-difference representation of this term is employed 
in the marching region, so analytic forms of these derivatives are only required in the 
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nose region. The following algorithm is employed. The shock values for p and v in 
the nose region are fit in a least-squares sense with the following conic 


Al e + B{(S- So) + CS (S - So) + 1 - = 0 

OO 


(4.3.6) 


where S is either p s or v s . This equation is constrained to pass through the stagnation 
value with a slope of zero. The shock gradients are supplied by the derivative of this 
expression: 

dS 2A£ + B(S-So) 

= £ 1 ( 4 -3.7) 

4 C(2S-S 0 ) + Bt-± 

*~>0 

On the stagnation line, the equilibrium shock jump conditions of Section 4.3 sim- 
plify to become 

1 


% = 


P*Q 


Ps o — Poo 1 




hsQ ^OG Hh 


Js 0 — lst)(PsQy 


K-i) 

~ o (P 5 0 ’ ^ 3 0 ) 


p* 0 


7 V T 

/v arr J Ay 


(4.3.8) 

(4.3.9) 

(4.3.10) 

(4.3.11) 

(4.3.12) 


Equations Eqs. (4.3.8) through (4.3.12) are solved using the iterative procedure of 
Section 4.3. The stagnation line value of the strearnwise derivative of u a is given in 
Section 3.3. The quantities 


(-*0 I u 3 d( 


[ u 3 d£ J 


are required for Eq. (2.4.13). Recall that the gradients are given by 


dS _ ±JHS- l Sq) 

d( C(2S-So) + B{-±- 

*So 

where S is either p a or v s . Also note that 

€- , ° l “s d£ J k Sq ( r , d£ J 


(4.3.13) 


(4.3.14) 
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Since r 3 


Further, 


so that 


4.4 


£ and S — > S 0 as £ — > 0, the above equation gives 


AS 

since AS/ A£ — * dS/d£ — + 0 as £ — * 0, 




Surface Quantities 


(4.3.15) 


(4.3.16) 


(4.3.17) 


The relations given in Section 3.4 for perfect gas flows can be employed for equi- 
librium flows as well, provided Pr eq is used in Eq. (3.4.8). 
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5 Nonequilibrium Analysis 

In nonequilibrium flows, the characteristic time scales of the chemical reactions 
are of the same order as the time scales of the mean flow. Thus, species conservation 
equations must be solved in order to determine the local composition of the fluid. In 
addition, it is advantageous to rewrite the energy equation in terms of T. The rate of 
production of each species appears in the energy equation when it is written in this 
form. These changes to the governing equations presented in Chapter 2 are detailed 
in the sections that follow. Further, discussions of the thermodynamic and transport 
properties, state equation, boundary conditions, and surface quantities are presented 
in this chapter. 


5.1 Viscous Shock Layer Equations 

The continuity and streamwise momentum equations and Maslen’s expression, 
presented in Chapter 2, are still valid for these flows. However, it is advantageous 
to recast the energy equation in terms of T (see Appendix G for details), since the 
chemical rate equations are functions of T: 



u_dT dT\ 
in ds +v !m) 


u Op dp 

h | 8 s dn 


/- 

■ dT’ 


dn 


+ k 


dr 

dn 


(j_dhi 1 (du 

\/ti dn ^ hj dn J ^ ^ 


u dh\ 
hi dn 


(5.1.1) 






i=l 


dT 

dn 


J^hiW, 

t — 1 


The subscript i is the species index for the seven-species air model used here. Table 5.1 
provides the key for this indexing, along with the species molecular weights and perfect 
gas values of specific heat. 

Each species present in the mixture is governed by a species conservation equation 
(see Appendix G): 


P 


( u dc{ 

~hllh 



- W, 


, IdJi (i dh, i dMl 
‘ \ dn +J ' t*. dn + h 3 dn )j 


where the binary diffusion mass flux is 


(5.1.2) 


J,= 



Lt i2 


dci 

dn 


(5.1.3) 
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Table 5.1. Species Data 


i 

Species 

Mi 

[gm/gm-mole] 

[cal/gm-mole-K] 

i 

n 2 

28 

7.00 

2 

0 2 

32 

7.00 

3 

N 

14 

4.97 

4 

0 

16 

5.44 

5 

NO 

30 

7.00 

6 

7 

N0 + 

e~ 

30 

0.0005486 

7.00 


and Le i2 is the binary Lewis number. Binary diffusion assumes that the fluid con- 
sists primarily of heavy particles diffusing to light particles, and vice versa. Previous 
investigations [8, 13] have used this approach successfully (assuming Lt\i = 1.4 ev- 
erywhere, for example). In keeping with the approximate nature of this approach, a 
constant Lewis number of 1.4 is used here. The species mass fractions are 


Ci - 


El 

P 


and by definition satisfy the relation 


(5.1.4) 


Nm 

S> = i 


(5.1.5) 


t=l 


In the normalized shock-oriented coordinate system (see Chapter 2), Eqs. (5.1.1) 
and (5.1.2) can be written as 

energy: 

PapCp 


l u 3 u 

dT r hl -\dn b dT' 

v dT\ 

\h r 

d£ n b d£ dr} n 

n b dr} n J 


U S U 

X 


dp rfn — 1 dll b dp 
d( Ps n b d£ dr} n 


+ 


v dp 
n b dg n 


l 8 

' dT 

\%i 

[ d Vn J 


+ 


dT l 1 dh\ ]_dhA 2 ( &* « gjtiV) 

' dr}, , \hi dr}n + hzdTjn) + sfi hi dr) n ) J 


c 2 N, 


dT 


N, 




1 = 1 


(5.1.6) 
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species conservation: 


U 3 U \dc x //„ - 1 dn b Oa 1 dc t 


hi [ d£ n b i)q n J n b di) n 


w + i! / JL [ 7. jM , f dc ‘ ( 1 1 ^3 

\%i [ ‘%J l di}n\h x dq n h 3 dr) n 


where 


k Le Vi dti _ Ji dc { 
C p , n h dt) n n b d rj n 


Ji — -7r~Lei2 — -jr Leu 

C Pf Pr 


(5.1.9) 


In preparation for numerically integrating the governing equations, the production 
terms are written in terms of 7’ and c, in the energy and species conservation equation, 
respectively. These representations are given in the next section, along with the 
resulting governing equations. 


5.2 Species Rates of Production 

In addition to their presence in the species conservation equations, the species 
rates of production («;,•) appear in the energy equation when it is written in the 
above form. These source terms are functions of both T and c, (see Appendix H). 
As such, they should be expressed in terms of c, for the species continuity equations, 
and in terms of T for the energy equation [58, 59]. 

For the species continuity equation, the production term can be linearized [58, 59] 
to yield 

— = w t - aw] (5.2.1) 

where 


= 

1! 

< -3 


+ r 

.> r:,,] 

p 

(5.2.2) 

h’i = 

Mi * 

^rej r = 1 

'KrKr 

+ r 

C, 

j 

p 

(5.2.3) 

r,+ = { 

° 5^ 

'-i 

i 

0 

T , for 
, for 

ft, 

ft. 

- <*,> 
~ °i, r 

> 0 

< 0 

(5.2.4) 

r- / 

a/, iT - 

r , for 

o, >r 

- fi i<T : 

> 0 

(5.2.5) 

'"“l 

0 

, for 

<*i,r 

- ft> : 

< o 
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and the remaining quantities are given in Appendix H. The subscript p denotes 
evaluation at the previous iteration. The quantities o, ir and are the stoichiometric 
coefficients for reactants and products, respectively. 

For the energy equation, the production term can be written as a truncated Tay- 
lor’s series [58, 59]: 


W{ Wi 

P P 


+ 


d_ 


[T - T r ] 


(5.2.6) 


where again the subscript p denotes evaluation at the previous iteration. Thus, the 
production term in the energy equation can be expressed as 


N, 


= w l{ + Tw 2i 


i=l 


where 


and 


Further, 


x 


*-{'&•»(?)}. 

d ( wA _ \ 

dT\p)~ T ~ ^ 1,r a, ’ r ^ 

\b,, + ^ - a.) R„ - (b,„ + ^ - a) A.r 


(5.2.7) 


(5.2.8) 


(5.2.9) 


(5.2.10) 


Now make use of these expressions in the energy and species continuity equations. 
For clarity, the entire set of governing equations for nonequilibrium flows is repeated 
below. 
continuity: 


d , _ - 1 , V n 1 dnf) d , _ _ , . ps d / - l l \ n 

- {p.pu,uhz) - (MJ - (pvhlhz) = 0 


(5.2.11) 


£ - momentum: 




dii _ du s Ti n — 1 dn b du 


U *dt + 6 di Us iib d( drin J n b [di In ' hr drf n 


du u dhr 

+ 


1 

'dp 

Tjn ~ 1 drib 

dp 

e 2 

( d L | 

f du 

u dhA 

h\u 9 


rib d£ 

dr) n . 


[dr]n r 

\9rtn 

^1 J m 


+ P 


1 2 dhr 1 dhA ( du u 

, h i dr/n h 3 dr) n ) \di] n A dr) n ) J 


(5.2.12) 
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energy: 


PaPCp 


(u a u\d'T T) n — 1 dn b dT 


Pf \ h i n b d( drj n J n b drj n 


u a u dp Tj n — 1 dn b dp v dp 
h i pi n b d( dip, n b drj n 


t 2 \ d r dT] , , dT ( i din i dh 


n b { d'h L OVn J fyn \hidp n h :i dtp, J 

c 2 N ‘ dT 

+ ^ J,Cpt d^ t ~ " , ‘ 1 “ Tw ' 2 ' 


du u dh \ \ 2 
dr) n hidijn) 


species conservation: 


u 3 u \dci i] n - 1 dn b dc, ] v dc{ 


hi [di n b di di) n I n b drj n 


Maslen’s equation: 


, K*r a u s v, sin 

P = Pa + — — (V - 1) 


KaVs 1 /„ 2 


cos r, 


(v 2 - 1) 


•o -ii , e f d j. dci ^ dc, f 1 dhi 1 dh 3 \] 

Vi-CiWi Ji~r— + J7,— — + 1 -^? 5.2.14 

J n 2 b \di] n dip, drp L \hi di) n h 3 di) n J ) v ’ 


(5.2.15) 


= n l P* u l J o Ph (Vn 1) dp n ] + n b [p s u s r a £" pudp n (5.2.16) 

Eqs. (5.2.12), (5.2.14) and (5.2.14) are cast in the standard form for parabolic 
partial differential equations to yield the following coefficients (see Appendix G): 
i - momentum: 


A 0 = ~ r (l 


1 nj dip, fl \hidp n h 3 dp„ 


„ [ u u s u Tj n - 1 dn b 

~ Psp h -r — 

n b hi n b di 


£ 2 / 1 dh x \ r dp 


A'} 9 ( , 


1 din 1 dh 3 


n l \h\ drpj [dtp, f \/i, dip, h 3 dip, j J 

A _ 1 17 1 jgl Vn - 1 dp f 1 dn b 


_ ( _du s v dhi 

p3p \ di n b di) n 


hi Uu.df J n b dtp, \u s d£ 


A 4 = p a u 3 
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energy: 


Ay — 


nt 


A\ — — - 


n 


b L 


aJb . ( l dh\ l dh 3 \ ^ ~ 

wC k {ir,wCT 3 wJ +n ^ JiC, \ 


dr) 

~PsP 


A 3 — — 


U S 'U 

TT 


•y « s u ?/„ — 1 </n(, 

rib It i n >> d£ 

A 2 = U'2, 

dp 1) n — 1 d/l;, dp 

»lfc tfyn 


c, 


*7 


(5.2.18) 


v dp 
n b dr) n 


2 u 2 ( du u dhi \ 

efi n 2 b \dT) n hi dt) n J 


+ tii. 


species conservation: 


a P™ 

A x = p s u s — 
h\ 


Aq — — ~ J7i 
*6 


Ax = 


7^ 


9 


d*7n 


+ 3 ' {I'M + VM) 


+ P»P 


v u a ii i) n — 1 drib 

rib hi n b d£ 


A 2 = -p,pu>; 

Aj = p s pu>° 


i 


(5.2.19) 


A 4 p s ti. 


pu 

*7 


All iY, species conservation equations need not be integrated to determine the 
concentrations of each species. It is possible to use only N 3 — 1 of the species continuity 
equations, in conjunction with Eq. (5.1.5). Ref. [16] goes one step further by locally 
conserving the elemental composition of the fluid. For air, this gives two relations, if 
binary diffusion assumed: 


cjv 2 + cat + 
co 2 + c O + 


M 


N 


Mno 

Mo 


{cnO + C NO + ) = Cyv. 


200 


(' c NO + Cwo+) = C0 2c 


(5.2.20) 


(5.2.21) 


Mno 

so that only N a — 2 species continuity equations must be integrated. Numerically 
integrating fewer equations should reduce CPU requirements, so this latter approach 
is employed here. 


28 


5. NONEQUILIBRIUM ANALYSIS 


5.3 State Equation 


As mentioned in Chapter 4, K is variable for reacting flows since M of the 
mixture is not constant. When finite-rate chemical reactions are considered, the 
flowfield can be modeled as a mixture of thermally perfect gases. For such a mixture, 
the equation of state is 

Hu w 

p = p- 


M 


where 


Hu = 


3L 

Qoo 


The mixture molecular weight can be written as 


M = 


1 


N a 


Ci 




(5.3.1) 


(5.3.2) 


5.4 Thermodynamic and Transport Properties 

5.4.1 Thermodynamic Properties 

The thermodynamic properties are required for each species that is present in the 
fluid. Expressions from Ref. [60] are used to define Cp i and /(, for the temperature 
range of 300 A’ < T* < 30000 A' (see Appendix I). The thermodynamic properties for 
the gas mixture are determined in terms of the individual species properties through 
the relations 



N, 

h - 

1 = 1 

(5.4.1) 

and 

Ns 



L P j = 

;=:i 

(5.4.2) 


5.4.2 Transport Properties 

Ref. [60] also supplies the species transport properties (//, and A:,) with curve fits 
for the temperature range of 1000A < T* < 30000A (see Appendix J). While the 
thermodynamic properties of the fluid are linear combinations of the individual species 
thermodynamic properties, the transport properties are defined by more complicated 
expressions. These expressions are referred to as mixing laws. Ref. [61] presents a 
thorough review of mixing laws, and points out several approximate formulations of 
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them. In the present research, two of these simplified relations are employed. The 
method of Armaly and Sutton [62] is used to define the mixture viscosity: 


where 


and 


N, 


yid — — + 


N, ,.2 

t=l 


XiXj MiMj 


V, ~ fUj (Mi + MjY [3 An Mi 


5 + M, 


(5.4.3) 


(5.4.4) 


/^u “ r 




VFl + VJd 13 


12 


(5.4.5) 


The mole fraction for species i is 


X { c, 


M 

Mi 


(5.4.6) 


Further, Aij = 1.25 and 

f 0.78 , for interactions of c 
— 1 0.15 , for interactions of e 


atoms or molecules with each other 
atoms or molecules with ions 


(5.4.7) 

For the mixture conductivity, the method of Mason and Saxena [63J is employed: 

k ‘ (5.4.8) 


M a 

1=1 


Ms 


' + 


J = 1 


where 


1.065 


I + 


<b = 

' 3 2^2 


[PjMl 

n m7 




(5.4.9) 


Ref. [61] shows that values from the curve fits approach Sutherland‘s values as 
T* 1000 /v. Thus, for temperatures less than 1000A' , Sutherland’s viscosity law 

/i = 1.4584 x 10 — — 11Q 33 
and Sutherland’s law for thermal conductivity [32] 

^★1.5 


[gm/cm-s] 


(5.4.10) 


k* = 5.9776 x 10“ 6 - 


cal/cm-s-K\ 


T * + 194.4 

may be used to define the mixture viscosity and conductivity, respectively. 


(5.4.11) 
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5.5 Boundary Conditions 

5.5.1 Wall Values 

As in the case of perlect gas and equilibrium (lows, the no-slip condition is applied 
at the wall so that 

= v w = 0 

Further, the wall temperature, T w {£), must be specified. Additional boundary con- 
ditions must be supplied for the species continuity equations. For a noncatalytic 
surface, the boundary condition is 


dci 

dr) n 


\W 


- 0 


If equilibrium catalytic wall conditions are specified, then 


(5.5.1) c 


^tw Cj e<l (l w ) ( 5 . 5 . 2 ) 

In this research, the wall temperature is low enough to produce freestream values for 
these concentrations: 


c 0 , w = 0.23456 c N . iw = 0.76544 (5.5.3) 

C0 W = C N W = <'N0 W ~ t'Afo+ u, = v e - w - o 


5.5.2 Shock Properties 


for nonequilibrium calculations, flow through the shock wave is assumed to be 
chemically frozen. The governing equations for the shock jump conditions axe (see 
Appendix F for details): 


v 3 = 


sinf, 

Ps 


(5.5.4) 


Vs = Poo A Sin 

* 2 i '' 

sin 1 




h s — h ^ “j- 


k / 8 > 


C (hoo h 0 2 (^5) ~h CN 2 ^ h yv 2 ( T s ) 

Cvs = ^ 00 C f0a (T.) + cs 2 JJ PNi (T,) 

_ Ps 
Ps - JO 

u T a 


M 


(5.5.5) 

(5.5.6) 

(5.5.7) 

(5.5.8) 

(5.5.9) 
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where the species enthalpies of Eq. (5.5.7), as well as the species specific heats in 
Eq. (5.5.8), are obtained from the curve fits of lief. [60]. 

By substituting Eq. (5.5.9) into Eq. (5.5.5), the following quadratic for p„ is ob- 
tained (with T a as the only unknown): 

p 2 s - ( Poo + sin 2 r.) Ps + ~T 7 ~~T s sin 2 Y s = 0 (5.5.10) 

' •A'r air 


Substitute Eq. (5.5.9) into Eq. (5.5.6) to obtain 


hs I 




•W air Ps _ J 


(5.5.11) 


Differentiating Eq. (5. 


n u 


Sill 


2 r, 


9ps _ May. 

2 Ps - (pco + sin 2 r s ) 


Differentiating Eq. (5.5.11) with respect to T s gives 

Ohs (KuVT.tl T s 0p 3 

0T 3 ~ * m p.[p. PlOTs 

Combine Eqs. (5.5.12) and (5.5.13), and rearrange to get 

n u 

Ohs 


sin 2 T s 


OTs 


PsPs - sin 2 I\, 


(5.5.12) 


(5.5.13) 


(5.5.14) 


As with equilibrium flows, these equations must be solved in an iterative fashion. 
To speed convergence, Newton’s method is employed: 



Ts = T Sp + A T a 

(5.5.15) 

where 

^ | r ^ 

1 

II 

3 

(5.5.16) 

The function F is defined to be the difference in enthalpy values given by Eqs. (5.5.6) 
and (5.5.7), respectively: 


h [Is] - ha Eq (5 5 G) ^ s Eq. (5.5.7) 

(5.5.17) 

so that 

F\T.) = -C„ 

(5.5.18) 


The solution procedure is as follows: 
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1) Begin with 1\ from the previous shock evaluation, 

2) Use Eq. (5.5,10) to calculate p 9 . 

3) Determine a new value of p s from Eq. (5.5.9). 

4) Find h, from Eq. (5.5.6). 

5) Evaluate li, from Eq. (5.5.7) and C p , from Eq. (5.5.8). 

6) Calculate F and F' from Eqs. (5.5.17) and (5.5.18), respectively, 
to determine AT, from Eq. (5.5.16). 

7) If | AT, | is greater than the prescribed tolerance, update T, using 
Eq. (5.5.15) and return to (2). 

8) After convergence, compute v, from Eq. (5.5.4). 

I he derivative of u, is given in Section 3.3. The derivatives of p, and v, with 
respect to £ are (see Appendix F) 


dp. 


K . COS r. 


{v, - p, sin r.) 


a 


— 9 


Pa 


T 

— -f- 




C, 


PS J 


(u a - sin Fs)| 


■A4 a it 

7Z U 



and 


respectively. 


do, 

~di 


2k, cos V, 



1 

Ps 


r , 1 d Ps 

K, COS r s + — — 

V, d( 


(5.5.19) 


(5.5.20) 


5.6 Stagnation Line 

As mentioned in Chapter 2, a limiting form of the governing equations can be 
obtained for use on the stagnation line. A limiting form for the energy equation 
cast in T, as well as the species conservation equations, may be found in the same 
manner. These ordinary differential equations are given in Section 5.6.1. In addition, 
the stagnation line jump conditions and their gradients are presented in Section 5.6.2. 


5.6.1 Viscous Shock Layer Equations 

continuity: 

2K So n bi) p SQ piih l - P, 0 ~ (pvh’i) = 0 (5.6.1) 

£ - momentum: 


P^P 


( u 2 du, 

h i d( 


v 



du 

Oljn 


u dhi 
hi dlln 


+ — lim 
A, f-o 



dp 


y n — 1 drib dp 
n b dt, dp n 
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c 2 r & 

n bQ \ dt] n 


( du u dhi 

dr) n hi drj n )\ 


+ 3 


p dhi ( du u dhi \ 1 

'ti dr} n \dr)„ hidrjnjj 


energy: 


... or d P 


(5.6.2) 

(5.6.3) 


f 0 

' dT' 

\ di] n 

[ 


, dT 2 dhi] dT r . , 

+ / + £ - ni ° K + ru ”. 


species conservation: 

dci 

- A * = 

Maslen’s equation: 


w® — Ciibj 


. i!_/A 

«6 0 \ dif n 


Ji 


dci 


drjn 


’ d*?n /ii d?7„ 


v = Ps 0 ~y (»? 2 - 1 ) 


with 


i r / ,,n 

-r/ = n bQ 2 U 0 ^, 0 2 j o pu (f/„ - 1 ) d?/ r 


+ "6 0 K) K *0 


/‘ 7?n 

, / 

io J 


(5.6.4) 


(5.6.5) 


(5.6.6) 


Eqs. (5.6.2), (5.6.3), and (5.6.4) are cast in the standard form for parabolic partial 
differential equations to yield the following coefficients (see Appendix G): 

£ - momentum: 


A o — — ■ 


»*fcn 


A\ — — 


dp p dh\ 


A, = 


jL 

«6 0 2 \hl dl) n 

1 

As — -j— 

hi 


n bo 2 [dg n h\ drjn 


dp . 0 P dhi 


P’ 0 P 


«fcn 


[di) n ~ hi dfj n J 
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hi \ d£ rib 0 di] nt 

i—o | u, J dil n t~0Yu,dl JJ 


(5.6.7) 


energy: 


A. = 0 


A o — — ■ 


«6 n 


Ai = 


dk ^ k dhi 


,2 /V 3 


”6 0 2 L^„ 

A 2 = W 2; 


+ - p°oP„, Ci 


«6n 


*7 


(5.6.8) 


34 


5. NONEQUILIBRIUM ANALYSIS 


, i> dp 
A 3 = -X h 'U>1 
«6 0 OT) n 

A 4 = 0 


species conservation: 


»*n 


A, = 


<9 


(iTi) 4- 2 


^lo — — 

Ji dh\ 

I V " V ' “/il <9r/ n 

^2 = -Ps 0 pw] 

A'3 = P, 0 PW° 

A L, = 0 


+ V 


»6n 


(5.6.9) 


5.6.2 Shock Properties 


On the stagnation line, the shock jump conditions of Section 5.5.2 simplify to 
become 

1 


with 


% = 


P* o 


7\/ J - (Poo + l)p S0 + 


— Loo "t" 2 ^ 1 


-fe.“ T = 0 

A<«, J ° 


7; 


l2\ 


*0 


. air P^q J 

^ s o c 0-2 OO Lq 2 (t, 0 ) + c/v 2oo /ijv 2 ^7 jq^ 

^PaO — C °2 oo^P0 2 (^o) C ^2oo ^PyV 2 ( T *o) 

P’0 


P° 0 


7^-ti rp 

M atr 30 


Ru 


dh 3 

dT, 


M, 


(5.6.10) 

(5.6.11) 

(5.6.12) 

(5.6.13) 

(5.6.14) 

(5.6.15) 


r (5.6.16) 

P* 0 P> o ~ 1 

Equations Eqs. (4.3.8) through (4.3.12) are solved using the iterative procedure of 
Section 5.5.2. 
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Next, the limiting forms of the shock derivatives defined in Section 5.5.2 must be 
determined. The stagnation line value of the derivative of u a is given in Section 3.3. 
The limiting forms of Eqs. (5.5.19) and (5.5.20) are 


Km{±£ 

«-° 1 u s d£ 


1 V^- 2 [^ + c L lK-i) 


'so L^O o. 

M~r /T, 1 \ 

1 O T ^ 




(5.6.17) 


I u a d£ ) 


% / p * 0 L v *o l u ’ d Z 


respectively. 


(5.6.18) 


5.7 Surface Quantities 

The relation for skin friction which was defined in Section 3.4 is still valid for 
nonequilibrium flows. However, for chemically reacting flow, the effect of diffusion, in 
addition to that of conduction, must be included in the heat flux calculations. The 
heat transfer to the wall is 

C = k* + A (5.7.1) 

where q* d is the energy flux due to diffusion in the normal direction. This term can 
be expressed as (see Appendix G) 


N, q N , 

d = -e'VnZ'tTA = D«* 

»'=1 UU «=1 


Written in nondimensional form (see Appendix A), 


— -T . dc, 

~ t — e 2 _^ p r ^ ^ 


(5.7.2) 


(5.7.3) 


From Section 3.4, the heat flux due to conduction is 


2 dT t 1 dT 

q c = -e fr— = k- 

OU 71b 07) n 


Combining these components, 


f)T Na 

q w = -t 2 k—-^h t Ji 

i dn J, 


(5.7.4) 
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or 


dT ^pLen dc, 

"” = - £ l^ + g-pT^j 

For a noncatalytic surface, 

dci 


n b L rJfj n fr[ dq n J 


di/n 


0 


and tlie heat transfer rate reduces to the familiar form 

qw = -e 2 


9 

dr’ 

_ t 2 

, dT 

t 2 

k dn 

«, r| fr 

k rj 

<tyn\ 


(5.7.6) 

(5.7.7) 

(5.7.8) 
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6 Turbulence 

For turbulent flow, the governing equations are modified [41] using an eddy vis- 
cosity model. In various efforts [33, 41, 42, 11, 43, 44] the Cebeci-Smith [45, 46, 47] or 
Baldwin-Lomax [48] eddy viscosity models have been mated with VSL solvers. Such 
algebraic (or zero-equation) models are less complicated than more exact approaches, 
such as the Johnson-King [49] and two-equation [50] models, and as a result are more 
computationally efficient (although theoretically less accurate). The equations in this 
chapter may be applied to perfect gas or equilibrium flows. 


6.1 Viscous Shock Layer Equations 


The VSL equations of Chapter 2 are modified [41] using an eddy viscosity model 
so that they are applicable to turbulent flows. Simply put, the laminar transport 
coefficients are supplemented with turbulence terms. Thus, replace 


/i with fi + /q 

where /q is the eddy viscosity, and replace 

k with k + k t 


where k t is the eddy thermal conductivity. 
By definition, 

Pr = 

JL = 

Pr 


so that 


Cpfi 


k 

a 


The turbulent Prandtl number can be defined as 

C> 


Pr t = 


k t 


so that . 

k-t _ kt 

T^ l ~'c v 

Thus, the thermal conductivity may be rewritten as 


JL + ^L z= -L + J!± = -?-(\ + Pr \ 

Cp C p Pr Pn Pr\ fiPrt) 

In this study, a value of 0.90 is used for the turbulent Prandtl number [44]. 
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Summarizing, the modifications to the governing equations are 


P 


(l + C) 


_P_ 

Pr 


~( 

Pr V 


, Pr\ 

1 + £ Pn) 


where 


e + _ ti 


p 


These modifications affect only the energy and streamwise momentum equations. 
Referring back to Eqs. (2.1.2) and (2.1.3), these equations become 
s - momentum: 


( u du 

du 

uv dhi i 

\ l dp 2 

f d_ 


du 

u dhi 

\Vr7h 

+ 

<2 

4- 

hi dn ) 

J ~ 

' ^ hi ds 

( dn 

(l + C+J 

' dn 

^ hi dn 


+ 


2 dhi 1 dhfi 


hi On h 3 On 


l . . \ du u dhi 1 ) 

"( 1 + ' (6L1 > 


energy: 


( u dh dh\ u dp dp 

P TIT + v ^~ ~ v ~ = 


k\ ds 

\Ji\jL (i + e +fL\2!L] + JL( 

\ dn Pr\ Pr t ) dn Pr\ 


dn 


h i ds dn 




+/i(l +t + ) 


'duY 

d - 2 


u dhi du 
hi dn dn 


r\dh( 1 dhi 1 
\) dn \ Ax dn + h 3 

( u dhi 

hilfa 


dh 3 


dn 


( 6 . 1 . 2 ) 


Note: The described modifications are not performed on the higher order 

terms. 

Looking in the transformed system, Eqs. (2.1.11) and (2.1.12) become 
£ - momentum: 


PsP 


U 

Jh 


du _ du, r/„ — 1 dip, du 
' d£ + d( U " n b d£ drf n 


du u dhi 

+ 


n b [drf n hi dr) n 


1 

4 

'dp 

Tfn — 1 dn b dp 


e 2 f d 


IS 

hiu. 

k 

n b d£ drf n _ 


n b \ drf n 

ftl drf n \ 


+ 

/ 2 dhi 1 dh 3 

\hidtfn h-sdrfn 

) 

P (l + f+ ] 

1 du u dhi 

drfn ^ hi drfn 



(6.1.3) 


P.P 


(u,u 

dh ij n — 1 dn b dh 

v dh \ 


n b d£ drf n 

n b drf n ) 


energy: 
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U a U 

fit 


dp >jn — 1 drib dp 

aT 


+ 


v dp 


lib dlj n I n b dl} n 


(6.1.4) 



it 

( x Pr> 

i dh 1 

L l 

/ . Pr\ 

dh 

( 1 dhi 

1 dh A 

{ dt) n 

r 

Pr 

( 1 + £ 7v 

1 

C 

R* 

+ K- 

( 1 + e rFj 

1 drjn 

V K drjn 

h 3 dr) n ) 




0 + 


h x f ( du \ 2 u dfh du 1 ( u dhiV 1 

' \ \dp n ) “hi dr) n drjn J ^ \hi dp n ) j 


As a result, Eqs. (2.1.14) and (2.1.15) become 
£ - momentum: 

/t„ = -^(i+ e + ) 

- 2 r o 


Aj — ~ 


ni 


dr) n 


W ' + t+ )} + , ‘( 1 +c 


2 dhi 1 dh 3 

+ 


-pep 


n b 


+ 


[hidrjn h 3 dijn t 
u s ii ijn — 1 drib 


p dh\ 
hi dr) n 


hi rib d£ J 


(6.1.5) 


1 dl,, 


ni 


III dl) n 


An — 


dp 

dljn + f \ h 


1 


hi 


1 dhi 1 dhA] psp_ l _dus _ v dhi 

h[dij n + hi drjn )\ hi \ d£ n b di) ni 

1 dp 1 ijn — l dp ( 1 drib 
is d£ J n b dljn \ u a d£ 


A,t = p s u 


f)U 

% 


energy: 


e 2 pf , Pr\ 

Ao = ~ri~FrV +t ~P7j 




= -~( 
nl V 


2 ' ,Pr\ 

1 + £+ d 


-psP 


+ 


_dij 

u a ii ij n — 1 drib 


iib 


lb lib d£ 


L (JL\ 1 dhi 1 dh 3 \ 

<n \Pr) Pr\hidr) n h 3 dq n ) 

£i "-f (> 

rdij n \ 


n 2 Prdrjr 


l+e 


PrJ 


Ax = 


ll 3 ii 


—t 


2< 

nl 


u dhi 
hi drj n 


A 2 = 0 

dp ?/„ — 1 drib djp 

d£ n b d£ drjn 




( 6 . 1 . 6 ) 


+ 


-2 


t; dp 
rib drj n 

u dhi du 
hi di) n drjn 


A 4 = p,u s 


pu 

Ti 


Note: For e + — 0 , the laminar expressions of Eqs. (2.1.11), (2.1.12), 

(2.1. If), and (2.1.15) are recovered from Eqs. (6.1.3) through (6.1.6), 
respectively. 
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6.2 Stagnation Line 

If fully turbulent flow is to be considered, the limiting forms of the governing equa- 
tions are needed in order to calculate stagnation line solution. As £ — ► 0, Eqs. (6.1.3) 
and (6.1.4) reduce to 
£ - momentum: 


P>oP 


( u 2 du a v 
hi d£ n bQ 


du it dh\ 

+ 


+ 


t 

energy: 


dr)n hi dt] n J 
Mj_f (i . +\ jidhi' 

0 2 \ dt )n L ' 9, ln ^ hi 9 rj n 


T~ lini j — 

hi o ( u s 


+ 


3 dhi 

hi drj n 


dp jj„ - 1 dn b dp 

d( 

P 


n b d( dr) n 

du 
drj n 


} 




( 6 . 2 . 1 ) 

u dhi 


jLIJL ' 

n bQ \drj n 


JL (\ + c + —\ — 
Pr V Pr t ) dp n 


hi dg n J j 

dr)n ' ' dr, n ~ (6,22) 

2 dhi p / + Pr\ dh | 

hi di) n Pr V Pr t J drj n J 


dh dp 
P’i)P v — + = 


It follows that Eqs. (6.1.5) and (6.1.6) become 
/ - momentum: 


Aq — — - 


™b n 


(> + f+ ) 


Ai = 


A 2 = 


n b 0 
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d 

b 

1 dhi 


dp , A p dhi] t ^ / du, _ i> dM 
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+ 
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A,= 


hi 


lim { — ^7} - — — lim j — 

[c-o ( u, d( J n fc() drj n i-o ( u, d{ J 


0 


energy: 


'“-~tp?( l+e+ K;) 


A\ 


=--( 
«6n V 


1 -f t 


+ f* / /* \ . 3 dhi p 


-ps u pv - 


Pr t ) [ dfj n \ Pr) + h x drj n fr\ 

-f^( l+£+ ^ 

n 6 0 Prdijn \ 

a 2 - 0 


Pr t J 


(6.2.4) 


A3 = v 


dp 

drjn 


A 4 = 0 
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6.3 Turbulence Model 

A turbulence model must be employed to define e + , and the Cebeci-Smith model 
is used here. It is a two-layer eddy viscosity formulation whose inner layer value, ef , 
is based on Prandtl’s mixing length concept [51]. The outer layer value, , is given 
by the Clauser-KlebanofF [52, 53] expression. The inner layer expression is used from 

the wall outward until 4 > 4 - 

Note: The Cebeci-Smith model was developed for a body-normal system, 

however it is being applied here in a shock-normal environment. 

In the laminar to turbulent transition region, the composite eddy viscosity is modified 
using the Dhawan and Narashima [54] method. 


6.3.1 Inner Layer 

The inner layer eddy viscosity is based on Prandtl’s mixing length theory: 


4 = 


pi 2 l 
t 2 fi nb 


du u dh\ 


\drj n h x di] v 

The mixing length for the inner layer according to Van Driest [51] is 

/ = K v n b i]n 1 - ezp (-n + /A + )] 
where the von Karman constant is 

1< V = 0.4 

and the normal coordinate parameter is 


n + = 


nbi]„p 

f&W 

du 



pn b 

dl In 

w- 


1/2 


In Ref. [44], the damping factor is defined as 

*-»(£r 

where the local shear stress is 

T ( 1 + 1 ) |^ n h x dr) n 
The subscript u, indicates that quantity is evaluated at the wall. 


du u 3h\ 


(6.3.1) 


(6.3.2) 


(6.3.3) 


(6.3.4) 


(6.3.5) 
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6.3.2 Outer Layer 


The outer layer eddy viscosity is approximated by the Clauser-Klebanoff [52, 53] 
expression: 

(6.3.6) 


c + 
" O 


t; = 0.0168- 


e 2 /i 


The displacement thickness is 


h = n b 



(6.3.7) 


where 6 is the value of 
factor [55] is 


i) n at the boundary layer edge. 



Klebanoff’s intermittency 


(6.3.8) 


6.3.3 Boundary Layer Edge 

By definition, at the boundary layer edge, 


H oo 

where the total enthalpy (H) is 


H = h + 


u 2 -f V 2 


Numerically, this can be approximated as the first grid point where 

II 


H 0 


> 0.995 


This approach can lead to oscillations at the boundary-layer edge. Reference [56] 
proposes a different criteria, which is based on another characteristic of the boundary 
layer edge: 

d / II \ 

dijii 

Numerically, this is approximated as 

d ( H \ , 

d Vn \hJ ~ °’ 5 

This is the criteria which is employed in the present work. 
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6.3.4 Transition 

Transition from laminar to turbulent flow is modeled through the definition of a 
streamwise transition intermittency factor [54]: 

= 1 — exp(— Q.412^ 2 ) (6.3.9) 


where 


t = 2.96 


(| - fr) 

[to(X - 1)1 


(6.3.10) 


Before transition begins, the flow is laminar and 7 i,( = 0. As £ * 00, it can be seen 

that 7;^ -> 1. In practical applications, this value of unity is essentially achieved at 
£ = 5. Thus, transition begins at £ = £0 ( as prescribed by the user) and continues 
until £ = 5. Within the transition region, the local value of eddy viscosity is modified 


by 


t + = t + 7.,« 


(6.3.11) 


The value of X is dependent on the flow conditions and body geometry under 
consideration [57]. It is usually determined empirically. In the results presented 
later, a value of A = 2 is chosen for the purposes of comparison. 
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7 Shock Layer Geometry 

This chapter discusses the relations which define the shock and body shapes, 
and the resultant shock layer thickness. Sections 7.1 and 7.2 discuss the analytic 
expressions used to describe the body and shock shapes, respectively. Relations for 
the shock layer thickness and its streamwise gradient are presented in Section 7.3. The 
final section contrasts three coordinate transformations which have been investigated 
as a part of this research. 


7.1 Body Geometry 

7.1.1 Forebody 

For many configurations of interest, the body shape can be defined by an axisym- 
metric conic equation: 

r l =2(z b - n hQ ) - b b (z t - n t() ) 2 (7.1.1) 

where n bQ is the stagnation standoff distance and the nondimensional nose radius of 
curvature has a value of unity. This is the equation of an axisymmetric conic body 
whose nose is at ( z b = n b ,r b = 0). 

The character of the conic is dictated by b b : 

<0 , hyperboloid 

, = 0 , paraboloid 

b > 0 , ellipsoid 

. = 1 , sphere 

For a hyperboloid, the value of the cone angle, F fr , that the hyperboloid asymptotically 
approaches is specified so that 

b b = — tan 2 Ft (7.1.2) 

For b b > 0, the value of b b is a required input. By setting b b > 0, Eq. (7.1.1) can be 
used to define the forebody of an elliptic-cone, where b b = 1 is the special case of a 
sphere-cone. 

For the blunted body described by Eq. (7.1.1), the slope is 

(iv b 1 b b iS n bQ ) 

dz b 


r b 


(7.1.3) 
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7.1.2 Stagnation Line 

From Eq. (2.1.4), the body radius is 

r b = h 3b = r a - n b cos T* = r a |l - n b 

On the stagnation line this becomes 

v b — v s 1 v a h 

so that Eq. (7.1.3) can be rewritten as 

dr b 1 — b b — « 6 0 ) _ 1 ^ 

dz b r s h\, r a h\ 


cos r. i 


(7.1.4) 


(7.1.5) 


7.1.3 Ellipsoid- Cone Juncture 

At the juncture between the ellipsoidal lorebody and cone afterbody, both the 
position and slope of the body are continuous. Thus, the coordinates of this juncture 
may be determined by equating Eqs. (7.1.3) and (7.1.11). Solving for z b = Zj un ct 
(after some manipulation) gives 

Ft | (7.1.6) 

6 + sin 2 T b J 


Now the corresponding value for r Junc t can be found from Eq. (7.1.1). 



I'junct ~ ~ ( Z ]unct ~ 71 b 0 ) h ( Zj unct ™6 0 ) 

(7.1.7) 

For the special 

case of a sphere-cone ( b b = 1): 



Zjunct — ^ 

(7.1.8) 

and 

Rjunct COS r 6 

(7.1.9) 


For z b < z junch the body is defined by the equations of Section 7.1.1. After this 
juncture, the body shape is described by the equations of Section 7.1.4. 


Zjunct — 11 b 0 1 i V 


sin 


yj b b COS 2 T, 


7.1.4 Afterbody 

The cone afterbody is given by 

r b — i'j unc t T (~6 Zjunct) tan T^ (7.1.10) 

where rj unct and z ]unct are the coordinates of the ellipsoid-cone juncture, defined in 
Section 7.1.3. The slope of the cone is 

dr b 
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7.2 Shock Shape 

7.2.1 Subsonic- Transonic Region 

The body shape in the nose region is described by a conic equation (see Sec- 
tion 7.1.1). Van Dyke and Gordon [64] suggest that a conic body generates a conic 
shock shape, so the subsonic-transonic region is described by a conic equation: 


r- 2 — f 


2C a z, - 


(7.2.1) 


where C s is the value of the shock radius of curvature on the stagnation line and b s 
determines the conic’s character. 

Since the shock shape is defined analytically, closed-form expressions for the quan- 
tities listed below may be determined, (see Appendix K for details). The shock angle 
can be found from 

8 ir»r 3 = (c 3 -b 3 z a ) < ^ 


The curvature is given by 


cos 


Ac q ^ 


— (sin 2 T, -f 6, cos 2 F,j 


Further, the gradient of the curvature is 


i I k s cos 2 r, 

~jr — 3/Cj sin T, (6, — 1) 

d( r. 


(7.2.2) 


(7.2.3) 


(7.2.4) 


Referring to Eq. (2.2.4), the streamwise pressure derivative can be rewritten as 


dp 


dp, _ T] - J 


^ | 

4 ^ sin F, (u, + 3 it, cos 2 T, ( 6 , — 1) + /c,r,j 


dv. 


. ' ( 2 v , k , sin 2 r s cos r, (6, - 1) + (l + ) 

4 ( \ cos 1 , / 

In addition, from Eq. (2.2.7), the streamwise derivative of n becomes 


sin T,~ - v , k , cos T* 
dt. 


(7.2.5) 

) 


On 




= (*, - 1)-^ COS F, (6,-1) 


(7.2.6) 


On the stagnation line, Eq. (7.2.4) has the following limiting form: 

1 d/c, ) 


liiTi , „ 

o 1 cos T, d £ 


3s: s 2 (6, - 1) 


(7.2.7) 
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Thus, Eq. (2.4.12) becomes 


Km =Hin +^(^-1) 

C -0 1 u, d(, I €-o [ U s d( ) 


i] 2 — 1 


2 ~ V « K « (*-“2) + gm|— - 5 ^- 


and Eq. (2.4.16) becomes 




(7.2.8) 


(7.2.9) 


7.2.2 Supersonic Region 

Since the supersonic region is hyperbolic-parabolic in nature, a marching proce- 
dure is employed. The shock shape segment between the previous (k- 1) and current 
( lc ) stations is described by a truncated Taylor’s series (see Appendix L for details): 



r s — ro + T\ (Az*) + (Az 3 ) + ^ (A^,) 

(7.2.10) 

where 

A y — r — 7 

L-l*, s — ~ s 


The shock angle 

can be found from 



tan r, = — = r, + f 2 (Az a ) + — (Az,) 2 

CIZ g « 

(7.2.11) 

The curvature is 

given by 

3 r d 2 r, 
K s = - COS 1 5 

(7.2.12) 

where 

(P'r 

/: = >• 2+h(Az 3 ) 

dz 2 

(7.2.13) 


In this marching region, a simple backward difference representation of the pres- 
sure derivative is employed. Thus, at a given point k on station j , 


dp _ Pk, j - Pk-i,j 

8£ “ A <‘ 


where 


A & = & - £k - 1 
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7.3 Shock Layer 

With each iteration in the solution of the governing equations, a calculated value 
of the shock layer thickness may be determined from the continuity equation. In 
addition, the analytic expressions for the shock and body shape provide the means of 
calculating a geometric value of the thickness. If the equation used for the shock shape 
is reasonable, then these two values should match. This observation serves as the 
criteria for shock-shape convergence at a given station (see Chapter 8). Expressions 
for the calculated and geometric thickness values are given below. In addition, an 
analytic expression for the gradient of the shock layer thickness (a term which appears 
in the governing equations) is presented. 


7.3.1 Calculated Thickness 


The continuity equation can be used to determine the shock layer thickness based 
on the conservation of mass. This calculated thickness (see Appendix D) is defined 
by the following quadratic for n b : 


n l cos E, jf pu fa 7„ - 1) dr ] n j + n b J/> 8 u 4 j\, pudr ) n j - = 0 

Dividing by Eq. (7.3.1) can be written as 

2 . u 9 COS F s f l ] f u /I 1 

n 6 — 2 y o pd (>7n - 1) drj n + n b |2 Ps~ J q pudg n J -1 = 0 

As £ — > 0, Eq. (7.3.2) becomes (see Appendix E) 

n b J o pd (Vn - 1) dj] n +n b jf -pudj ] n J -1 = 0 

This expression is used to define the stagnation standoff distance, 7 i& 0 . 


(7.3.1) 


(7.3.2) 


(7.3.3) 


7.3.2 Geometric Thickness 

Through geometric considerations, 


r b = r, - n bg cos T 3 (7.3.4) 

where 

n b g = \[{z, - 2fc) 2 + (r, - r t ) 2 (7.3.5) 

is the geometric shock layer thickness (measured along a shock normal). 
Differentiating with respect to f gives 
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Similarly, 

~b ~ ~s + n b (J s i n 

Differentiating with respect to £ gives 

dm 


dz b . • r 

— ^cosl. + ^s.nf. 


(7.3.7) 


(7.3.8) 


However, 


(7.3.9) 


dv b dv b d^ b 

d£ dzb d£ 

where dr b /dz b is the local slope of the body geometry. Combine Eqs. (7.3.6), (7.3.8), 
and (7.3.9) to get 

An 


dn b 

~d£ 


. „ dr b 

Sin 1 a — cos I a 

dz 6 


, dr b . 

cos I J + - — sin 1 a 

d?b 


(7.3.10) 


which is valid for any body geometry. 
A shock normal has a slope of 


m 


'dr, 
, dz s 


-l 


Thus, the equation of the shock normal line passing through the shock point (z b ,r b ) 
and intersecting the body at point (^,7*6) is 

r b = rnz b + {r s - 1112 ,) (7.3.11) 


Forebody 

Combining Eqs. (7.1.1) and (7.3.11) gives 

( rn 2 + b b ) zi + 2 (r, - mz s ) rn — b b n b 0 — I z b + (r, — mz s ) + b b n bo + 2n bo = 0 
V ’ 1 (7.3.12) 

This quadratic can be solved for z b , and then r b is found from Eq. (7.1.1). With dr b jdz b 
defined by Eq. (7.1.3), the gradient of the shock layer thickness can be determined 

from Eq. (7.3.10). 

Stagnation Line 

The limit of Eq. (7.3.10) as £ -> 0 is 


1 dn b 


In 


1 dn r 

1 — COS 1 4 

dz b 


'"olcosl'a d£ / COsTa 

dz b 


(7.3.13) 
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Now combine litj. (7.1.5) with this result to get 


lim 


dn h 


cos F s d( 







so that 


lim 
£— "0 


1 drib 
cos r a d£ 



1 


(7.3.14) 


(7.3.15) 


Afterbody 

Combining Eqs. (7.1.10) and (7.3.11) gives 

— rjunct ~ ~ ( r * ~ mz ») 
m — tan 1\ 


(7.3.16) 


Substitute into Eq. (7.1 . 1 0) to get . Utilizing Eq. (7.1.11), the shock layer thickness 
gradient is determined from Eq. (7.3.10). 


7.4 Spacing Across Layer 

At the beginning of this investigation, a considerable effort went into the devel- 
opment of an inviscid nonequilibrium algorithm. In this early work, the governing 
equations were transformed to a streamline coordinate system (Figure 7.1). In this 
transformed space, the solution was advanced by marching along streamlines. At each 
new station, the shock streamline is added to the grid. Thus, if a constant number 
ol points is desired, an automated procedure for dropping streamlines is required. 

In Ref. [29], y (where // = ^/^s) was found to be a better transformation coor- 
dinate since it allowed more control of the spacing than the streamline method. In 
addition, this transformation gives a constant number of points across the layer for 
each station. Since viscous effects were to be included eventually in the nonequilib- 
rium algorithm, a switch from the streamline approach seemed inevitable. Unfortu- 
nately, this decision meant abandoning the explicit approach under development by 
this author. 

In Ref. [29], the solution is advanced along lines of constant y (Figure 7.2). The 
spacing is still based on mass flow through the layer, rather than on physical distance 
across the layer. While this transformation to y is an improvement over it still 
yields a grid-point distribution which is less than optimal in the near-wall region of 
the shock layer. In fact, the points cluster near the shock, rather than near the wall. 
As will be shown in Chapter 9, the shock layer profiles from the method of Ref. [29] 
exhibit oscillations as a result of this spacing. 
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Figure 7.3. Equally-spaced lines of constant i] n . 


in older to better control the spacing across the layer in the present technique, a 
transformation to normalized coordinates as described in Chapter 2 is used. Recall 
the definition of the following computational variables: 




Vn 



n b 


This normal spacing based on Tj n ( 1' igure 7.3) is analogous to the spacing used in 
Ref. [9]. Although Figures 7.2 and 7.3 show equally-spaced lines, these spacing func- 
tions lend themselves to clustering. However, since the physical distance to which 77 
corresponds varies in the streamwise direction, determining a distribution which will 
work well for the duration of a given solution can be difficult. On the other hand, 
since there is no streamwise variation in the normalized distance that Tj n represents, 
finding a suitable distribution is more straightforward in the present approach. 

The following distribution function [9] may be used to define the spacing across 
the shock layer: 


h~ ft) 

T /nj — JO + 7 r ~JT 

1 +cxp[ri, lp [r/ UM - (j - l)jj 


(7.4.1) 


where 


fo = 


Vn 


Jac 


J - * 


Vn 


[ Jac 


3a— 1 



(7.4.2) 
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(7.4.3) 
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In addition, j, is the total number of grid points, r/„ trt is the point through which the 
exponential fit passes, 7 lnj ac is th e geometric stretching factor, t)n ini is tbe intercept on 
the axis which determines the slope of the linear distribution near the shock, and 
is the slope of the fit at ij ntrk . This function clusters points according to a geometric 
progression near the wall, a linear distribution near the shock, and an exponential 
fit for the interior of the shock layer. If 7/„ in< < 0 is specified, a simple geometric 
progression is used so that 


dljuj — 7 lnj ac df)rij_ i 

Vn } - + drj nj 


and 


(L)nt ~ 


Vn 


l Jac 


- 1 


Ti i 3 ^ — 1 
* n fac 


(7.4.4) 

(7.4.5) 

(7.4.6) 


This relation clusters points near the wall, but not at the shock. 

For the perfect gas and equilibrium solutions presented in Chapter 9, j 3 = 101, 
while for the nonequilibrium solutions, j 3 = 51. The factors corresponding to these 
values of j 3 are 


is = 51 

’hlrk = 30 
Vn fac = 1-17 

Vn iKt = 0. 

'K lr = 0.13 


Js = 101 

r}n brk = 60 

V-fac = h13 
1" ml = O' 

T N >lp = 0- 1 
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8 Method of Solution 

The procedure for solving the equations presented in the preceding chapters is out- 
lined here. Fiist, an introduction to the issue of generating the shock shape is given. 
Next, the quasi-linearization of the streamwise momentum equation is discussed. The 
procedure for numerically integrating the equations of standard parabolic form is pre- 
sented in Section 8.3. Ihe chapter concludes with a discussion on the advancement 
of the solution from one station to the next. 


8.1 Shock Shape 

The full VSL approach requires the user to input an initial shock shape, which is 
used to calculate a first iteration solution. The shock shape calculated in the first- 
iteration solution must then be smoothed and used as the input shape for a second 
iteiation. dhis user-controlled process is continued until the new calculated shock 
shape varies little from the input shape (usually 2 to 4 iterations). In the design 
environment, such requirements of the user are undesirable. 

Automating this process was addressed in Ref. [29] by incorporating a shock it- 
eration technique in the solution. Since the subsonic-transonic region is elliptic in 
nature, this portion of the flowfield must be solved in a global fashion. The deriva- 
tive of shock layer thickness is described in terms of a cubic equation involving the 
body radius where an iterative procedure is applied to determine its coefficients. This 
process typically requires ten to twenty iterations to obtain a converged shock shape. 

A marching technique is employed aft of the subsonic-transonic region, since the 
inviscid layer is supersonic. For this region, the shock shape at the current station is 
extrapolated from the previous station using a cubic equation for the shock radius as 
a function of axial location. Requirements that this expression match the position, 
slope, and curvature at the previous station leaves one free coefficient (essentially the 
shock curvature derivative). An iterative technique is employed to determine this 
value based on comparisons between the shock standoff distance obtained from the 
flowfield solution and the geometric distance between the shock and body curves. 
When these two values agree to within the prescribed tolerance, the solution is ad- 
vanced downstream to the next station. The iterative procedure used in Ref. [29] is 
sensitive, so this marching technique is often slow to converge. 

The current technique also generates its shock shape as part of the solution. The 
details are given in Section 8.4. Before discussing that topic, the linearization pro- 
cedure is given in Section 8.2. Fhen, an outline of the solution of the governing 
equations is presented in Section 8.3. 
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8.2 Streamwise Momentum Linearization 


The usual approach to handling the nonlinearities in the VSL equations is to eval- 
uate the A m coefficients (rn = 0, 1, 2, 3, 4) from the most recent information available. 
For example, the value of u obtained in the most recent solution of the streamwise 
momentum equation is used in loading the coefficients for solving the energy equation. 
This approach, which essentially treats all such quantities as knowns, still leaves the 
following nonlinear terms in the streamwise momentum equation: 

du _ 9 


du 

'Wn 




For the energy equation, the nonlinearities are handled through a simple lagging 
technique. However, in order to speed convergence, the streamwise momentum equa- 
tion is quasi-linearized [65]. This process involves expressing u as 


u D + A u 


( 8 . 2 . 1 ) 


where the subscript P denotes evaluation at the previous iteration and A u is the 
change in u from the previous iteration. Differentiating Eq. (8.2.1) with respect to 

' 7ngiveS - - ■ /a -\ 

m 


du du 


+ A 


dl) n dl] n 

and differentiating Eq. (8.2.1) with respect to £ gives 


( 8 . 2 . 2 ) 


du 

ae 


du 

ae 


+ A 


(8.2.3) 


Using these representations, 

du 

dr) n 
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so that if the higher order term is dropped, then 
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d Vn P d Vn 
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(8.2.4) 


Using Eqs. (8.2.1) and (8.2.2) this can be rewritten as 
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(8.2.5) 
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Similarly, 
and finally, 
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( 8 . 2 . 6 ) 


ir = 2 {wp ii} 


ul 


(8.2.7) 


1 he “bracketed” portions of these three expressions are the terms which would be 
present without linearization. Eqs. (8.2.5) through (8.2.7) are substituted back into 
the ^-momentum equation to obtain the linearized form of the equation. 

This linearization process only affects the A 2 and A 3 coefficients of the stream- 
wise momentum equation. Based on Eq. (6.1.5), the coefficients for the streamwise 
momentum equation become 
£ - momentum: 
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On the stagnation line, these become 
£ - momentum: 
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(8.2.9) 
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8.3 Equations of Standard Parabolic Form 

A finite-difference approach is used to solve the ^-momentum, energy, and species 
conservation equations (for nonequilibrium flow). As mentioned earlier, the first step 
is to write these equations in the standard parabolic form: 


d 2 W dW . . dW n 

Aq — — + A\— 1- AiW + A 3 + A 4 ^— — 0 


rivl 


rirjn 


dt 


(8.3.1) 


where W represents the dependent variables w, h, and c,'. Next, finite-diflference 
representations are substituted for the partial derivatives. A two-point backward 
difference is used for the partial derivatives with respect to £. For a fully implicit 
scheme, a three-point central-difference of values at the current station is used for the 
partial derivatives with respect to rj n . Thus, for an ai’bitrary point J at station k 


ri 2 W _ 2 [Wkj+i - (1 + /?) W kJ + fiWkj-i] 
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where 
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In a more general approach, the r/ n - derivatives are differenced about point ( k 1+0, j) 
where 

0 , fully explicit 

< 1 , Crank Nicholson 

1 , fully implicit 
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The derivatives now become 
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and 


where 


+ [WW, - (1 + 0)W„.,j -0)} 
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A?/n J+1 + P 2 AT) n . 
In addition, the quantity VE is evaluated at point (it - 1 + 0,j): 

w = w kJ e + w k - ld (i - ©) 

Making use of these expressions, Eq. (8.3.1) takes the form 

AjWkj-\ + Bj W k<J + CjWkj+i = Dj 
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On the stagnation line (since A 4 = 0 for the streamwise momentum, energy, and 
species conservation equations), the coefficients reduce to 
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The results presented in Chapter 9 were generated with the fully-implicit approach 

(0 = 1). 

Evaluating the coefficients of Eq. (8.3.8) at discrete points across the shock layer 
yields a tri-diagonal system of equations which may be solved using Thomas’ algo- 
rithm [66]. The general solution is 


W kJ = EjW ktJ+l + Fj 


(8.3.9) 
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where 
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The boundary condition at the shock is 
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On the body, the general boundary condition is 
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In finite difference form, this can be written as 
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this gives 


Ei = 1 


F\ =0 


(8.3.17) 
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Implementation of the algorithm involves first evaluating Ej and Fj from j = 2 up to 
J = Ja ~ I-, and then solving Lq. (8.3.9) for from j — j a — 1 down to j = 1. 

For a given station and shock shape, the governing equations are solved in the 
following order. The streamwise momentum and energy equations are solved for u 
and h, respectively. The equation ol state gives p. From the continuity equation, 
the shock layer thickness is determined. The pressure distribution is found directly 
from Maslen’s expression. The normal velocity distribution is extracted from the 
continuity equation after the standoff distance has been determined. 


8.4 Advancing the Solution 

The previous sections in this chapter explain the solution procedure of the govern- 
ing equations for a known shock shape at a given station. More precisely, holding the 
shock geometry (and the resulting jump conditions) constant, the fluid equations may 
be satisfied in an iterative fashion. From this solution, a calculated shock layer thick- 
ness for the current station is determined. This value is compared to the geometric 
shock standoff distance, with the difference between the two dictating the adjustment 
of the shock shape. After convergence of the shock shape for a given station, the so- 
lution is advanced to the next station, using the current profiles to initialize the new 
station. The discussion now turns to how the shock shape is actually determined. 

Maslen’s original method [18, 21] is an inverse technique. That is, a shock shape 
is prescribed and the resulting calculated body shape is compared with the desired 
body shape. The shock shape is then adjusted in an iterative fashion until the desired 
body shape is obtained. Riley and DeJarnette [67, 68, 69] automate this process for 
blunt bodies where separate approaches must be used for the nose region and the 
remainder of the flowfield. 

The current approach parallels that of Ref. [29] in that the shock shape is generated 
as part of the solution, rather than being an input required of the user. However, 
the approach of Riley and DeJarnette serves as a model for the technique employed 
here, since it provides a smooth shock shape within relatively few iterations. In the 
subsonic region, a conic equation is used to describe the shock radius as a function 
of axial position. The marching procedure again utilizes a cubic equation with one 
parameter to be determined iteratively. The details of the iterative procedures used 
to determine the shock shapes for the subsonic and marching regions are presented 
in Sections 8.4.1 and 8.4.2, respectively. 

8.4.1 Subsonic- Transonic Region 

Since the subsonic- transonic region is elliptic in nature, this portion of the flowfield 
must be solved in a global fashion. As discussed in Section 7.2.1, the shock shape for 
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the subsonic-transonic region is described by 


r] = 2 C s z a - b a z 


(8.4.1) 


where C 3 is the value of the shock radius of curvature on the stagnation line and b a 
determines the conic’s character. For a given body shape, three parameters ( b a , C a , 
and Ub 0 ) constrain the geometry of the shock-body tandem (and hence, the shock 
layer). The value of n bo is the shock layer thickness (standoff distance) calculated 
from the stagnation line solution (see Section 7.3.1). 

The quantities b a and C s are determined through a quasi-Newton iterative pro- 
cedure. With each variation of these two parameters, the flowfield is solved for the 
entire subsonic region. The fluid equations are solved at eleven discrete stations be- 
ginning on the stagnation line (giving n bQ ). The values for the calculated shock layer 
thickness at two stations near the end of the subsonic region are compared with the 
values dictated by the geometry (see Section 7.3). Shock iterations are continued 
until the calculated values of n b at these two stations match the geometric values. 
Typically only four to six iterations are required to converge the shock shape. 


8.4.2 Supersonic Region 

Since the supersonic region is hyperbolic-parabolic in nature, a marching proce- 
dure is employed. Recall from Section 7.2.2, the shock shape for the region between 
the previous and current stations is described by a truncated Taylor’s series (see 
Appendix L for details): 
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where the shock derivative is given by 
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(8.4.2) 


(8.4.3) 


Written in this form, the only unknown is the second derivative of r a with respect to 2 , 
(which is proportional to the shock curvature) at the current station. This parameter 
is determined through the iterative process described in the next paragraph. Note 
that Ref. [29] iterates about the derivative of the shock curvature, which is a more 
sensitive procedure. 

The position, slope, and curvature of the shock at the previous station are known. 
Thus, the curvature at the current station is the only unknown in Eqs. 8.4.2 and 8.4.3. 
As a first guess, the d 2 r a /dz 2 term at station k is extrapolated from its values at the 
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previous two stations. With the shock geometry (and corresponding jump conditions) 
constrained, the fluid equations are solved. Then the calculated and geometric values 
of rib are compared to determine an error (^ err ). After perturbing the guess for 
^ r */d* a (which alters the shock position and slope), the governing equations are 
solved to determine a new value for 6 err . Now the second derivative at station k may 
be updated using the secant method: 


1 

"T3 

<Pr, 

S„r J 

f <Pr. 

d?r a \ 
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dz\ 

dz> 

p ~ ^ctt |p | 

1 dzi ~ 

dz] 

J 


The subscript p denotes a value from the previous iteration. Through successive 
applications of the secant method (accompanied by a solution to the fluid equations), 
tferr -*• 0, with shock shape convergence typically obtained in three to four iterations 
per station. 
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9 Results and Discussion 

Results generated with the present algorithm are given in this chapter, which is 
divided into three sections: perfect gas, equilibrium, and nonequilibrium solutions. 
In an effort to evaluate the accuracy of this new approach, the comparisons outlined 
below are made. Heating rates are shown alongside experimental or flight data, where 
available. In addition, the perfect gas solutions are compared with the results of the 
VSL [9] and Grantz [29] approaches. Comparisons with the VSL solver of Ref. [9] are 
also made for equilibrium flows. The nonequilibrium results are compared with those 
of Ref. [14]. Specific areas to be addressed are: 

1) surface properties, including examination of the stagnation region; 

2) shock shapes; 

3) shock layer profiles of dependent variables; 

4) run times. 

The algorithms employed here are fully-implicit (see Section 8.3). For the perfect 
gas and equilibrium results, 101 points across the layer are used by each of the methods 
(see Section 7.4). The nonequilibrium solutions are based on 51 points across the layer 
for the VSL and present method. The results are plotted in nondimensional form (see 
Appendix A), with the dependent variable on a logarithmic scale in most cases. Note 
that in this chapter, the superscript * is omitted on the dimensional quantities. 
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9.1 Perfect Gas 

Two cases are presented in this section. The first case is a 5° sphere-cone with 
a length of approximately 140 ft. Two solutions with different nose radii are calcu- 
lated in Ref. [70] in a study of bluntness effects. Comparisons with those results are 
made here. The geometry of the second case is a 15° sphere-cone with a length of 
approximately 20 in. Again the solutions are obtained for two different nose radii, 
and the results are compared with experimental data [71]. For each case, analysis of 
the second radius is limited to surface properties in the interest of brevity. 

Case 1 

For Case 1, Mach 15 flow over a 5° sphere-cone is calculated using the VSL [9], 
Grantz [29], and present approaches. The freestream conditions are p rXj = 0.018919 
psi and = 478° R, with a wall temperature of T w = 2259.69° R. Fully-laminar 
solutions are calculated, as well as solutions where the flow transitions to turbulence. 
The beginning of transition is an input for all three approaches. The VSL heating 
results are from Ref. [70]. ° 

First, the solution for R noae = 1.5 in is calculated for a body length of s b = 1125. 
The transition point is specified to be s b = 192. Figure 9.1 shows the body shape 
along with the shock shapes (fully-laminar) calculated with the three methods. The 
approximate (Grantz and present) techniques agree with the VSL results as the sharp 
cone solution is approached. However, from the close-up given in Figure 9.2, note that 
the two approximate approaches yield a thinner shock layer than the VSL algorithm 
in the pressure overexpansion region. 

In Figure 9.3, the computed heating results from the three solutions are presented 
for both laminar and turbulent flows. The laminar results of the VSL and present 
method are in excellent agreement overall, while those of Ref. [29] are at an approx- 
imately ten to fifteen percent higher level. Downstream of the nose region (for s b 
between approximately 10 and 100), both approximate techniques over-predict the 
VSL heating results. It is believed that Maslen’s pressure relation is the source of this 
deviation, as will be discussed in the next paragraph. The turbulent results of the 
VSL and present method are in excellent agreement overall, with those of Ref. [29] 
approximately twenty percent higher. The transition heating lesults (see Figure 9.3) 
of Giantz [29] reflect an error in the application of the transition model. 

Figure 9.4 shows the body pressure distributions from the three methods for the 
fully-laminar solution. The computed pressure distributions are in excellent agree- 
ment (within five percent) except in the pressure overexpansion/recompression region. 

In this region, the shock and body angles are quite different (see Figure 9.2) for slen- 
der cones, violating one of Maslen’s assumptions. Apparently, this poor pressure 
prediction feeds back into the shock shape which is generated as part of the solution, 
to give the shock shape deviation mentioned earlier. As a result, the VSL shock layer 
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Figure 9.3. Heat transfer comparison for 5° sphere-cone, R noae = 1.5 in. 



Figure 9.4. Body pressure comparison for 5° sphere-cone, R noae = 1.5 in. 
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thickness (and also the thickness gradient) in this region differs from the values given 
by the approximate methods. Since the thickness gradient plays a prominent role in 
the governing equations, this deviation is reflected in the shock layer profiles, as well 
as the surface properties. The effect of the differences in the overexpansion region di- 
minishes further downstream, as can be seen in the figures. A similar phenomena has 
been seen in the inviscid technique of Riley and OeJarnette [68], which also employs 
Maslen’s pressure relation. 

Figures 9.5 and 9.6 provide more detailed information on the stagnation region 
results. As the stagnation line is approached (s6 -+ 0), dq w ldst> -+ 0 and dp w /dsb — ► 0. 
Note that the heating results from the current method approach the stagnation line 
more smoothly than those of either Ref. [29] or the VSL method. This is a direct result 
of the limiting form of the governing equations used here (see Section 2.4). Note also 
that the results of Ref. [29] are approximately twenty percent higher than the other 
two solutions in this region, and do not approach a zero gradient at the stagnation 
line. In Figure 9.6, the stagnation region pressure distribution is reasonably smooth 
for all three methods. 

Skin friction results are presented in Figures 9.7 and 9.8. The method of Ref. [29] 
does not compute this parameter, hence it is excluded from this comparison. Both 
the laminar and turbulent results are in excellent agreement, (generally within five 
percent). This result is not surprising since Reynolds analogy relates the skin friction 
to the heat transfer. Therefore, good agreement for the skin friction is expected in 
those areas where the heat transfer is predicted well. As shown in Figure 9.8, the 
present method has a smoother distribution in the stagnation region. 

Now consider the stagnation line profiles. Figure 9.9 shows the pressure predic- 
tions to be within 0.5 percent of one another. In Figure 9.10, note that the normal 
velocity results of the VSL and present method are virtually identical, while the re- 
sults of Ref. [29] differ from them slightly. The approximate relation for v used in 
Ref. [29] is the source of this difference. The tangential velocity is zero on the stagna- 
tion line, but a finite non-zero value for its normalized form (u = u/u 3 ) is calculated 
from the limiting form of the streamwise momentum equation (see Section 2.4). Pro- 
files of this normalized velocity are shown in Figure 9.11. The approximate results 
are in close agreement with each other, while the VSL results are different. This dif- 
ference may be due to the VSL limiting form of the streamwise momentum equation 
(see Section 2.4). Figures 9.12 and 9.13 show that the enthalpy and density profiles 
of the three methods are virtually identical on the stagnation line. The density mir- 
rors” the accuracy of the other results (since it is calculated from the state equation), 
so in the remainder of this chapter it will be shown only sparingly. 

Pressure profiles near the pressure minimum (sj, ~ 40) and at the end of the body 
(sj, ~ 1125) are presented in Figures 9.14 and 9.15, respectively. The approximate 
results are virtually identical. Agreement with the VSL solution is excellent as the 
sharp cone limit is reached. Even in the pressure minimum region, where the deviation 
is largest, the Maslen pressures are generally within ten to fifteen percent of the VSL 
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Figure 9.5. Heat transfer comparison for 5° sphere-cone, 
Rnose = 1.5 in (stagnation region). 
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Figure 9.6. Body pressure comparison for 5° sphere-cone, 
Rnoa e =1.5 in (stagnation region). 
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Figure 9.9. Pressure profile comparison for 5° sphere-cone, 
Rnoae = 1.5 in (stagnation line). 
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Figure 9.10. Normal velocity profile comparison for 5° sphere-cone, 
Rnose =1.5 in (stagnation line). 



9.1. PERFECT GAS 


p„ = 2.72 psf 
T„ r 478 °R 
= 15 

T w = 2260 °R 


VSL 

Grantz 

Present 


o.oo l— 
o.oo 


Figure 9.11. Tangential velocity profile comparison for 5° sphere-cone 
R noae = 15 in (stagnation line). 
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Figure 9.13. Density profile comparison for 5° sphere-cone, 
R-nose — 1 ■ 5 in (stagnation line). 


results. 

For purposes of comparison, the velocities calculated in the two approximate 
methods (solved in a shock-normal system) are resolved into body-normal compo- 
nents. The velocity component normal to the body is shown in Figures 9.16 and 9.17. 
The limitations of the approximate expression for v which is used in the method of 
Ref. [29] are readily apparent. This component is small relative to u and h except in 
the near-wall region. As a result, its poor prediction has minimal effect on the overall 
shock layer solution outside of the boundary layer. In the pressure minimum region, 
the other two techniques agree well near the wall. The last station exhibits excellent 
agreement between the VSL and present approaches. 

The body-tangential velocities are shown in Figures 9.18 and 9.19, and the three 
approaches are in good agreement. The enthalpy profiles of the present and VSL 
approaches, shown in Figure 9.20 and 9.21, also agree well with one another. At 
the end station, the results of Ref. [29] differ in the boundary layer. This difference 
is reflected in the heating rates, and is due to the inaccuracies of the approximate 
relation for the normal component of velocity. Also at the end station, note the 
severe oscillations around the boundary-layer edge which axe present in the results 
of Grantz [29]. Figure 9.19 shows wiggles in the tangential velocity profiles as well. 
This is the result of not having enough points in this region. A grid adjustment is in 
order, but since the Ref. [29] normal spacing is based on the stream function rather 
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Pressure profile comparison for 5° sphere-cone, 
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Pressure profile comparison for 5° sphere-cone, 
Rnose = 1.5 in (s b ^ 1125). 
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Figure 9.16. Normal velocity profile comparison for 5° sphere-cone, 

Rnose — 1.5 171 % 40). 



V 

Figure 9.17. Normal velocity profile comparison for 5° sphere-cone, 
Rnose — 1.5 tTl ^ 1125 ). 
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Figure 9.20. Enthalpy profile comparison for 5° sphere-cone, 

— 1.5 1TI % 40). 



Figure 9.21. Enthalpy profile comparison for 5° sphere-cone, 
Rnoae — 1.5 ill 1125). 
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Table 9.1. Run-times 0 for 5° cone, R noae = 1.5 in. 



VSL 

Ref. [29] 

Present 

stations 

643 

333 

^ 209 

CPU time 

4232 

867 

408 

grid pts/sec 

15 

39 

52 

shock iterations 

3.0 

6.5 

2.3 

grid pts /sec/shock 

46 

252 

121 


a - Sun Sparcstation lH- 


than a physical distance, an adequate distribution is difficult to determine. These 
wiggles, which are documented in Ref. [27], were repeatedly observed early on in this 
research and were the motivation for switching from rj- to ^-spacing, which is based 
on distance across the layer. 

Since the nature of this research is to develop an approach for use in preliminary 
design, where computational speed is an important consideration, some comparisons 
in run time are presented. Table 9.1 shows the overall run times required to generate 
the turbulent solution for the 5° sphere-cone. The value shown for the VSL results 
is the total time required for three shock iterations. However, the CPU requirements 
to generate the initial and second shock shapes are not included, so this figure is a 
low estimate for the VSL method. It can be seen that the current algorithm yields 
a dramatic decrease in total CPU requirements over both the VSL and Ref. [29] 
techniques. Since varying numbers of streamwise stations are used by each of these 
approaches, the average number of grid points solved per second is also presented. As 
shown, three global iterations of the shock shape are required for the VSL solution. 
The approximate methods iteratively determine the shock shape as the solution is 
advanced from one station to the next. In this case, for a given station, the method 
of Grantz [29] requires an average of 6.5 iterations to converge the shock shape, while 
the present technique averages 2.3 iterations. The last entry in the table takes into 
account the number of shock shapes calculated, and gives the processing capabilities of 
each approach for a given shock shape. This illustrates that solving the approximate 
VSL equations of the present method is inherently 2.5 times as fast as solving the 
full VSL equations. The method of Ref. [29] is able to process many more grid 
points/sec/shock shape than the current approach and this can be attributed to its 
more approximate governing equations. 

The solution with a nose radius of R n0 ,e = 9 in and a body length of Sb end = 180 
will now be discussed. The transition point is an input whose value is s b = 72. 
Because of this larger nose radius, bluntness effects are present for the length of the 
body. In the interest of brevity, the comparisons are limited to surface properties. 
The results for this body are similar to the results for the first 200 nose radii of 
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Figure 9.22. Heat transfer comparison for 5° sphere-cone, R nose — 9 in. 

the first body. Figures 9.22 through 9.24 show that the heating rates, laminar body 
pressures, and skin friction distributions follow the same trends cited above. That is, 
the results of the approximate approaches are in generally good agreement with the 
VSL solution. Likewise, the shock-layer profiles (not shown) follow patterns similar 
to those of the R note =1.5 in solutions. 
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Figure 9.24. Skin friction comparison for 5° sphere-cone, R 
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Figure 9.25. Shock shape comparison for 15° sphere-cone, R nose = 1.1 in. 

Case 2 

For Case 2, Mach 10.6 flow over a 15° sphere-cone is calculated using the VSL [9], 
Giantz [29], and present approaches. The freestream conditions are = 0.01915 
psi and = SS^OTO*#, with a wall temperature of T w - 540 ° R. Heating rate 
comparisons are made with Cleary’s experimental data [71], 

First, the solution for R no , e = 1.1 in is calculated for a body length of s;, = 20. 

Figure 9.25 shows the shock shapes from the three methods. As with the 5° cone, the 
two approximate techniques yield a thinner shock layer in the pressure overexpansion 
region than the VSL algorithm. The flare at the end of the shock shape of Ref. [29] 
reflects an incorrect interpolation for the last computational station. 

In Figure 9.26, the heating results from the three solutions are presented and 
are seen to compare well (generally within fifteen percent) with the experimental 
data, k igure 9.27 provides more detailed information on the stagnation results. As 
before, the heating results from the current method approach the stagnation line more 
smoothly than those of either Ref. [29] or the VSL method. 

The computed pressure distributions (Figure 9.28) are in excellent agreement 
(within five percent) except near the sphere-cone juncture. This deviation is not 
as large as that seen for the 5° cone, where the difference in shock- and body-angles 
is greater. As mentioned earlier, Maslen’s relation works best when the difference 
between these two angles is small, so the above improvement is expected. The skin 
friction comparisons of Figure 9.29 show that the VSL and present method are again 
in good agreement (generally within fifteen percent). 
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The stagnation line results are very similar to those given in Figures 9.9 through 
9.13, so they are not presented here. As before, profiles near the pressure minimum 
(s 6 ^ 6) and at the end of the body (s b « 20) are presented. Figures 9.30 and 9.31 
show the pressure profiles at these two stations. Figure 9.30 shows the results from 
the two approximate techniques are similar and generally within five percent of the 
VSL solution. The agreement is excellent between the approximate results, which 
are virtually identical, and the VSL solution as the sharp cone limit is reached (Fig- 
ure 9.31). At the pressure minimum region, the Maslen pressures are generally within 
five to ten percent of the VSL results. 

The velocity components normal to the body are shown in Figures 9.32 and 9.33. 
The behavior noted for the 5° cone is present here, although the difference between 
the results of Ref. [29] and the others is not as large as before. The body-tangential 
velocity (Figures 9.34 and 9.35) and enthalpy profiles (Figure 9.36 and 9.37) also 
follow the pattern set in the previous case. 

The solution with a nose radius of R n ost — 0.375 in and a body length of $b end — 50 
is now discussed. Because the nose radius is smaller than above, the bluntness effects 
(which are visible for the length of the body above) are less pronounced at large 
values of s b for this case. In the interest of brevity, the comparisons are limited to 
surface properties. Figures 9.38 through 9.40 show that the heating rates, laminar 
body pressures, and skin friction distributions follow the same trends cited above: the 
results of the approximate approaches are in generally good agreement (within five 
to ten percent) with the VSL solution. Further, the shock-layer profiles (not shown) 
follow patterns similar to the solutions with R n0 sc — 11 
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Figure 9.30. Pressure profile comparison for 15° sphere-cone, 
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Figure 9.31. Pressure profile comparison for 15° sphere-cone, 
H, lost =1.1 in (s b « 20). 
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Figure 9.32. Normal velocity profile comparison for 15° sphere-cone, 

R nose = 1.1 in ( s b « 6.5). 



Figure 9.33. Normal velocity profile comparison for 15° sphere-cone, 

Rnose = 1-1 i« (*6*20). 
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u 

Figure 9.34. Tangential velocity profile comparison for 15° sphere-cone, 

Rnoae — 1.1 (^6 ~ 6 . 5 ). 



u 

Figure 9.35. Tangential velocity profile comparison for 15° sphere-cone, 

R-nose — 1.1 (tt (S(, ~ 20). 
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Figure 9.36. Enthalpy profile comparison for 15° sphere-cone, 
Rnose = 1.1 in (fift « 6.5). 



Figure 9.37. Enthalpy profile comparison for 15° sphere-cone, 

Rnose =1.1 in ( s b RS 20). 
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Figure 9.38. Heat transfer comparison for 15° sphere-cone, R noae — 0.375 in. 



Figure 9.39. Body pressure comparison for 15° sphere-cone, R nose 


= 0.375 in. 
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9.2 Equilibrium 

Two sets of results for equilibrium flow are presented in this section. Both cases 
are for fully-laminar flow. The discussion follows the format of the perfect gas section. 


Case 1 

As a first case, flow over a 5° sphere-cone (#„„« = 0.114 in) is calculated using the 
the VSL [9] and present approaches. Both approaches use Hansen’s [35] equilibrium 
air model. The freestream conditions are p<x> = 0.067917 psi and T<x, = 438° .ft, with 
a variable wall temperature input. Fully-laminar solutions are calculated for a body 
length of s beni = 1350. Heating comparisons are made with the Reentry F flight 
experiment [72]. The VSL heating results aie fiom Ref. [73]. 

Figure 9.41 shows the shock shapes from the two methods. The results are analo- 
gous to those of the perfect gas, 5° cone solution discussed earlier, although the shock 
layer is not as thick here. In Figure 9.42, the heating results of the VSL and present 
method are in very good agreement with each other (within ten percent), although 

they both underpredict the flight data. _ . 

The body pressure distributions are in excellent agreement (within five percent) 
except in the pressure overexpansion/recompression region, as shown in Figure 9.43. 
Figures 9.44 and 9.45 provide more detailed information on the stagnation region, 
and show the results approach the stagnation line in a manner similar to the perfect 
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Figure 9.41. Shock shape comparison for 5° sphere-cone, R noae = 0.114 in. 


ga !i It 6 SkU1 fnctlon resu,ts of the two me thods are presented in Figures 9.46 

an 9.47. The two solutions are in good agreement, generally within ten percent. 

., , ^f atl0n line results for and p are presented in Figures 9.48 

thiough 9.52, respectively. The two sets of results exhibit better agreement than 
is seen for the perfect gas case. Profiles near the pressure minimum (s b « 50) and at 
the end of the body (s b m 1350) are presented in Figures 9.53 and 9.54. The pres- 
sure profiles are within ten to fifteen percent of each other at the pressure minimum 
while they are nearly indistinguishable (less than five percent apart) at the end of the’ 
body. 1 he body-normal velocity components at the pressure minimum (Figure 9 55) 
are ,n agreement near the body, although further out in the layer there are noticeable 
differences. At the end station (Figure 9.56), the agreement is excellent. In fact 
the present results are smoother than those of the VSL. The body-tangential veloc- 
!ty (Figures 9.57 and 9.58) and enthalpy profiles (Figures 9.59 and 9.60) agree quite 
well (within approximately five percent of each other). Again note the oscillations 
in the VSL results for the end station. These wiggles appear in the present solution 
as well, and reflect a grid resolution problem at the boundary-layer edge. Since the 
normal distribution for each method is based on physical distance, this problem can 
be readily corrected by adjusting the spacing parameters (see Section 7.4). 

Table 9.2 shows the overall run times required to generate this solution The 
value shown for the VSL results is the total time required for three shock iterations. 
As with the perfect gas comparison, the approximate algorithm yields a dramatic 
decrease in total CPU requirements over the VSL technique. Further, the last entry 
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Figure 9.42. Heat transfer comparison for 5° sphere-cone, R noS e — 0.114 in. 
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Figure 9.43. Body pressure comparison for 5° sphere-cone, R n0 ae = 0.114 in. 
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Figure 9.44. Heat transfer comparison for 5° sphere-cone, 
llnose — 0.114 in (stagnation region). 



Figure 9.45. Body pressure comparison for 5° sphere-cone, 
= 0.114 in (stagnation region). 
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Figure 9.48. Pressure profile comparison for 5^ sphere-cone, 
Rnose 0. 114 in (stagnation line). 



Figure 9.49. Normal velocity profile comparison for 5° sphere-cone, 
R-nost = 0.114 in (stagnation line). 
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Figure 9.50. Tangential velocity profile comparison for 5° sphere-cone, 
Rnosc = 0.114 in (stagnation line). 



Figure 9.51. Enthalpy profile comparison for 5° sphere-cone, 
Rnose = 0.114 in (stagnation line). 
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Figure 9.52. Density profile comparison for 5° sphere-cone, 
Rnose = 0.114 in (stagnation line). 


m the table illustrates that solving the approximate VSL equations for equilibrium 
flow is inherently about sixty percent faster than solving the full VSL equations. 
This speed-up is not as large as the perfect gas value, which can be attributed to 
the complexity of the equilibrium air model relative to the perfect gas relations. 
This is because a large percentage of the total run time is dedicated to determining 
thermodynamic and transport properties, which are the same in both methods, so 
that the benefits of using the time-saving approximate equations are reduced. 
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Figure 9.55. Normal velocity profile comparison for 5° sphere-cone, 
f^Tiose = 0.114 iii % 50). 



Figure 9.56. Normal velocity profile comparison for 5° sphere-cone, 
Ftnosc — 0. 114 in ( s b « 1350). 
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Figure 9.57. Tangential velocity profile comparison for 5° sphere-cone, 

Rnose = 0.114 in (s b « 50). 



Figure 9.58. Tangential velocity profile comparison for 5° sphere-cone 

R nose = 0.114 in ls b » 1350). 
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Figure 9.59. Enthalpy profile comparison for 5° sphere-cone, 
Rnose = 0.114 in ( Sb 50). 



Figure 9.60. Enthalpy profile comparison for 5° sphere-cone, 
Rno.' = 0.114 in { Sb « 1350). 
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Table 9.2. Run-times" for 5° cone, R nose 0.114 ill. 



VSL 

Present 

stations 

667 

206 

CPU time 

4539 

714 

grid pt s/sec 

15 

29 

shock iterations 

3.0 

2.5 

grid pts/sec/shock 

45 

72 


a - Sun Sparcstation 1-h 


Case 2 

As a second case, comparisons between the VSL [9] and current method axe made 
for Mach 22 flow over a 35.5° hyperboloid. The nose radius is R no ,t - 3 - 4645 ^ ft 
and the fully-laminar solution is computed for a body length of s bend = 24.5. The 
freestream conditions are Poo = 0.00094 psi and T«, = 400° R, with a wall temperature 

of T w = 1998° fi. , , .. . 4 „ . 

Figure 9 61 shows that the shock shapes from the two methods he virtually atop 

one another. In Figure 9.62, the heating rate comparisons are very good ( within five 
to ten percent) for the entire length of the body. Likewise, the computed pressure 
distributions (Figure 9.63) are in excellent agreement (within two to five percent) 
everywhere. Referring to Figure 9.61, these excellent results can be attributed to the 
fact that with the large body angle, the shock- and body-angles are nearly equal for 
the entire length of the body (which more closely matches Maslen s assumptions). 
The skin friction comparisons of Figure 9.64 show that the VSL and present metho 

are again in good agreement (generally within ten percent). , 

Since the stagnation line results are similar to those given in Figures 9.48 through 
9 52 and the agreement between the solutions of the two methods is still excellent, 
they are not presented here. Again, profiles near the pressure minimum (a* « 15.5) 
and at the end of the body (s b » 25.5) are presented. Figures 9.65 and 9 66 show 
that the pressure profiles at these two stations are nearly indistinguishable (less than 
five percent apart). The only major differences in the body-normal velocity compo- 
nents (Figures 9.67 and 9.68) are near the shock, and thus have minimal impact on 
the surface properties. Further, the body-tangential velocity profiles (Figures 9.69 
and 9.70) and enthalpy profiles (Figures 9.71 and 9.72) are within approximately five 

percent of each other. 
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Figure 9.61. Shock shape comparison for 35.5° hyperboloid, R n0ie — 3.46457 ft. 



Figure 9.62. Heat transfer comparison for 35.5° hyperboloid, Rnoae = 3.46457 ft. 





Figure 9.66. Pressure profile comparison for 35.5° hyperboloid, 


R nose = 3.46457 ft {s b ^ 24.5). 
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Figure 9.67. Normal velocity profile comparison for 35.5° hyperboloid, 
R nose = 3.46457 ft (s b * 15.5). 
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Figure 9.68. Normal velocity profile comparison for 35.5° hyperboloid, 
Rnose = 3.46457 ft (s b « 24.5). 
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Figure 9.69. Tangential velocity profile comparison for 35.5° hyperboloid, 

R no , e = 3.46457 ft (s b ta 15.5). 
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Figure 9.70. Tangential velocity profile comparison for 35.5° hyperboloid, 

R noae = 3.46457 ft (s b « 24.5). 
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Figure 9.73, Shock shape comparison lor 6° sphere-cone, R nose = 1.5 in. 

9.3 Nonequilibrium 

Three sets of results for nonequilibrium How are presented in this section. The 
first is flow over a 5° sphere-cone with a noncatalytic wall. The next case is a longer 
20° cone, where two solutions are computed: noncatalytic wall and fully-catalytic 
wall. The VSL [14] and current method are employed for this analysis. The discussion 
follows the format of the previous sections, although additional profile comparisons are 
made. As a final case, the present method is used to calculate flow over a hyperboloid 
with both noncatalytic and fully-catalytic wall conditions. The results for this case 
are compared with Shuttle flight data. 


Case 1 

As the first case, comparisons between the VSL [14] and current method are made 
for Mach 25 flow over a 6° sphere-cone. The nose radius is R note = 1.5 in and the 
solution is computed for a body length of Hfnd = 5. The freestream conditions 
are p = 0.00794 psi and T = 486° R, with a wall temperature of T w = 2260 °R. 
Nonequilibrium flowfield calculations require additional boundary conditions for the 
species mass fractions, and a noncatalytic wall boundary condition (see Section 5.5.1) 
is specified here. The VSL heating results are from Kef. [70]. 

Figure 9.73 shows the shock shapes from the two methods. As with the perfect 
gas and equilibrium results, the computed shock shapes begin to deviate from one 
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another in the overexpansion region. Figures 9.74 and 9.75 show that the surface 
heating results are within fifteen percent of each other and approach the stagnation 
value quite smoothly. The pressure expansion is still underway at the end of this body, 
and Figure 9.76 clearly shows the differences between the VSL and present body pres- 
sures. Note that the pressures are approaching each other as the solution progresses. 
Figure 9.77 displays the same smooth behavior seen previously for stagnation region 
pressure distributions. The skin friction results of the two methods are presente m 
Figures 9.78 and 9.79. The two solutions are in good agreement (generally within ten 


The^ stagnation line results for p,v,ii,T, and p are presented in Figures 9.80 
through 9.84, respectively. The pressure predictions (Figure 9.80) are within approx- 
imately one percent of each other, while the normal velocity profiles (Figure 9.81 
agree to within five percent. The tangential velocity profiles (Figure 9.82) are not 
as close, but as explained earlier, this may be due to differences in the streamwise 
momentum equation on the stagnation line. Agreement between the temperature, 
as well as density, profiles (Figures 9.83 and 9.84) is excellent The species mass 
fraction profiles are presented in Figure 9.85. Note that for this figure the dependent 
variables are assigned to the vertical axis, rather than to the horizontal axis as is done 
for the other profiles. Again the agreement between the two methods is very good. 
Only trace amounts of ,V0+ are present for this case, so its mass fraction profile is 

not included in the figure. . 

Since the flow solution does not encompass the pressure minimum (see h igure 9. lb), 
only profiles at the end of the body (s b = 5) are presented. Figure 9.86 shows that 
the pressure profiles differ by twenty percent of each other at the end of the body. 
The body-normal velocity components (Figure 9.87) are in good agreement throug 
out the shock layer at this station. Notice that this velocity component is negative 
at this station, which is indicative of the streamlines spreading away from the body 
surface (see Figure 7.1). Although this behavior has not been shown previously, the 
same situation exists in this region for perfect gas and equilibrium flows about slen- 
der cones. The body-tangential velocity (Figure 9.88) profiles are in agreement near 
the body and shock, but deviate from one another (by as much as ten percent) in 
the interior of the shock layer. On the other hand, the temperature (Figure 9.89) 
and mass fraction profiles (Figure 9.90) agree quite well (within approximately five 

percent of each other). _. ft ni 

As a final comparison for this case, the electron density is considered. Figure 9.91 

displays the profiles for the stagnation line and the end station. Note the good 
agreement between the solutions. As the sharp-cone limit is reached, the normal 
spacing for the two methods are virtually identical. Thus, the oscillations near the 
shock for the end station (which are present in both solutions) are possibly due o 
grid resolution. Since they occur far away from the body, they have little effect on 

the surface properties. T , oira]ll . 

Table 9.3 shows the overall run times required to generate this solution. The value 
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Figure 9.77. Body pressure comparison for 6° sphere-cone, 
Rnose = 1.5 in (stagnation region). 
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Figure 9.78. Skin friction comparison for 6° sphere-cone, R„ 
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Figure 9.79. Skin friction comparison for 6° sphere-cone, 
Rnose = 1.5 in (stagnation region). 
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Figure 9.80. Pressure profile comparison for 6° sphere-cone, 
Rno,e = U5 m (stagnation line). 



Figure 9.81. Normal velocity profile comparison for 6° sphere-cone, 
Rnose = 1.5 ill (stagnation line). 
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Figure 9.82. Tangential velocity profile comparison for 6° sphere-cone, 
Rnose — 1 • 5 in (stagnation line). 



Figure 9.83. Temperature profile comparison for 6° sphere-cone, 
Rnose = 1 • 5 in (stagnation line). 
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Figure 9.84. Density profile comparison for 6° sphere-cone 
Rnose =15 in (stagnation line). 
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Figure 9.85. Mass fraction profile comparison for 6° sphere-cone 
Rnose = 1.5 in (stagnation line). 
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P 


Figure 9.86. Pressure profile comparison for 6° sphere-cone, 
Rnose = 1-5 in (sf, = 5). 



V 


Figure 9.87. Normal velocity profile comparison for 6° sphere-cone, 

Rnose — 1*5 m (^5 5). 



nonequilibrium 



Figure 9.88. Tangential velocity profile comparison for 6° sphere-cone, 

Rnose = 1-5 in (sb = 5). 



Figure 9.89. Temperature profile comparison for 6° sphere-cone, 

Rnose = 1.5 in ( S b = 5). ( 
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Tin 

Figure 9.90. Mass fraction profile comparison for 6° sphere-cone, 

R-nose — 1 Ll (.Sj, — 5). 



Figure 9.91. Comparison of electron concentration profiles for 6° sphere-cone, 

Rnose — 1.5 171. 
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Table 9.3. Run-times 0 for 6^ cone, R n0 se — 1-5 in 



VSL Present 

r stations 

46 

38 

CPU time 

759 

273 

grid pts/sec 

3 

k 7 

shock iterations 

3.u 

4.2 ^ 

grid pts/sec/shock 

9 

30 


a - Sun Sparcstation 1 + 


shown for the VSl results is the total time required or three shook .tout o^ A^ h 
the perfect gas comparison, the approximate algorithm yields a dramatic decrease m 
totaf CPU requirements relative to the VSL technique. Further, the last entry ml 
table shows that the the present method -.. process upP~x.ma^three t.mes as 
many grid points per second (for a given shock shape) as the VSL technique. 
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Figure 9.92. Shoc k shape comparison for 20" sphere-cone, R noae = 1.5 in. 

Case 2 

As a second case, calculations are performed for the Mach 25 flow over a 20° 
sphere-cone ( R nose = 1.5 in) with the same freest ream conditions as the previous case. 
Comparisons between the VSL [14] and current method are made for a body length of 
Sb 'nd = 250 - The freestream conditions are p /Xi = 0.00794 psi and = 4 86° R, with 
a wall temperature of T w = 22G0 °R. Solutions are generated for both a noncat’alytic 
and fully-catalytic wall condition (see Section 5.5.1). 

Figure 9.92 shows that the shock shapes (rom the two methods are indistinguish- 
able, once the flow has reached sharp cone conditions. Figure 9.93 shows that with 
this laiger cone angle, the difference between the shock shapes is still very small in 
the overexpansion region. From Figures 9.94 and 9.95, the approximate and VSL 
surface heating results are generally within five percent of each other for both wall 
conditions and approach the stagnation value quite smoothly. Note that in the nose 
region the effects of wall catalysis are large (more than fifty percent), while further 
downstream this difference diminishes (to approximately twenty percent). The results 
of the present method do appear to satisfy dq w jds\, = 0 at the stagnation line more 
closely than the VSL approach. The catalytic boundary condition has little effect on 
the pressure, so only the fully-catalytic results for the body pressure distribution are 
presented in Figure 9.96. In this figure, the values from the two methods are nearly 
identical for s h > 20. At the pressure minimum, the difference between results is 
approximately twenty percent. The behavior of the pressure in the stagnation region 
is identical to that of the previous case (Figure 9.77). The skin friction results of the 
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Figure 9.93. Shock shape comparison for 20° sphere-cone, 
Rnose = 1.5 in (overexpansion region). 


two methods for both catalytic conditions are presented in Figure 9.97. The two so- 
lutions are in good agreement for both wall conditions (generally within five percent), 
and show only a ten percent difference due to the wall catalysis. 

Since the freestream conditions are the same as the previous case, and the geome- 
try is still spherically-capped, the noncatalytic wall profiles for the stagnation line are 
identical to Figures 9.80 through 9.85. For the fully-catalytic w^l solution the pro- 
files of p, u, and u are nearly indiscernible from those of Figures 9.80 through 9^82, so 
they are not shown here. The temperature and mass fractions (Figures 9.98 and 9.99) 
are visibly different from their noncatalytic-wall counterparts, although agreement 
between the two methods is still within five to ten percent. The mass fraction profile 
of NO* is not included in Figure 9.99, since only trace amounts are present for this 

Again, profiles near the pressure minimum {s b « 3) and at the end of the body 
(s b J 250) are presented. Because the pressure and velocity profiles are not greatly 
affected by the wall catalytic condition, only the fully-catalytic wall results for the 
profiles of * v, and u are presented. Figures 9.100 and 9.101 show that at the pressure 
minimum the pressure profiles are within ten to fifteen percent of each other while 
the end station, the profiles are nearly indistinguishable. At the pressure rnimmum, 
the body-normal velocity components differ by twenty percent (Figure ■ )• S 

ure 9.103 shows these components to be in excellent agreement (within five percent) 
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Figure 9.94. Heat transfer comparison for 20° sphere-cone, R noae = 1.5 in. 



Figure 9.95. Heat transfer comparison for 20° sphere-cone, 
R-nose = 1.5 in (stagnation region). 
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Figure 9.96. Body pressure comparison for 20° sphere-cone, Rnose 



Figure 9.97. Skin friction comparison for 20° sphere-cone, R n ose - 
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T 

Figure 9.98. Temperature profile comparison for 20° sphere-cone, 
Rtwat: — 1.5 in (stagnation line). 





Figure 9.99. Mass fraction profile comparison for 20° sphere-cone, 
R n ose — 1.5 in (stagnation line). 
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Figure 9.100. Pressure profile comparison for 20° sphere-cone, 
Rnose = 1.5 in ( s b « 3). 
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Figure 9.101. Pressure profile comparison for 20° sphere-cone, 

R noae = 1-5 in {s b « 250). 
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Table 9.4. Run-tirnes u for 20° cone, R nose = 1.5 in. 



VSL 

Present 

stations 

587 

157 

CPU time 

7961 

525 

grid pts/sec 

4 

15 

shock iterations 

3.0 

2.1 

grid pts/sec/shock 

11 

33 


a - Sun Sparcstation 1 + 

at the end station. The body-tangential velocity profiles at the pressure minimum 
(Figure 9.104) are in agreement near the body and shock, but deviate from one an- 
other in the interior of the shock layer. Figure 9.105 shows that at the end station 
these profiles are virtually identical. Temperature profiles at the two stations given 
above are shown in Figures 9.106 and 9,107 for the noncatalytic wall, and in Fig- 
ures 9.108 and 9.109 for the fully-catalytic wall. At the pressure minimum, both 
sets of profiles are within approximately five percent of one another. The profiles at 
the end of the body are practically identical. 

Mass fraction profiles at these two stations are shown in Figures 9.110 and 9.111 
for the noncatalytic wall, and in Figures 9. 112 and 9.113 for the fully-catalytic wall. 
At the pressure minimum, both sets of profiles are within approximately five percent 
of one another. Differences in the profiles at the end of the body are indistinguishable, 
with the chemistry effects concentrated in the near-wall region. Only trace amounts 
of NO* are present, so those profiles are excluded from The figures. 

As a final comparison for this case, the electron density profiles for the three 
stations discussed above are considered. Figure 9.114 and 9.115 show generally good 
agreement between the methods for the noncatalytic and fully-catalytic wall solutions, 
respectively. The wiggles in the end station profiles at the boundary-layer edge are 
probably due to poor grid resolution there. 

Table 9.4 shows the overall run times required to generate this solution. The value 
shown for the VSL results is the total time required for three shock iterations. Trends 
from the previous case are observed again here. 
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Figure 9.102. Normal velocity profile comparison for 20° sphere-cone, 

Rnose = 1.5 in(.« 4 «3). 



Figure 9.103. Normal velocity profile comparison for 20 sphere cone, 

Rnose = 1.5 in (s b ~ 250). 
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Figure 9.104. Tangential velocity profile comparison for 20° sphere-cone, 

Rnose — 1.5 I'll .1). 



Figure 9.105. Tangential velocity profile comparison for 20° sphere-cone, 

R-nose — 1*5 ill % 250). 
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Figure 9.106. Temperature profile comparison for 20° sphere-cone, 

Rnose = 1 -5 *«• (> s 6 ~ 3 )- 
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Figure 9.107. Temperature profile comparison for 20° sphere-cone, 

Rnose = 1-5 M* ( S b ~ 250). 
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Figure 9.108. Temperature profile comparison for 20° sphere-cone, 

Rnose — 1-5 ITL * 3 ). 



Figure 9.109. Temperature profile comparison for 20° sphere-cone, 

H-nose — 1.5 tfl (s/> I-w 2o0). 
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Figure 9.110. Mass fraction profile comparison for 20° sphere-cone, 

R nose = 1.5 in ( s b » 3). 



Figure 9.111. Mass fraction profile comparison for 20° sphere-cone, 

Rnose = 1.5 in (Sfr « 250). 
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Figure 9.112. Mass fraction profile comparison for 20° sphere-cone, 

Rnose 1.5 IH (s(, ~ 3). 
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Figure 9.113. Mass fraction profile comparison for 20° sphere-cone, 

Rnoae — 1-5 111 ~ 250). 
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Figure 9.114. Comparison of electron concentration profiles for 20° sphere-cone, 

Rnose = 1-5 in. 



Figure 9.115. Comparison of electron concentration profiles for 20° sphere-cone, 

Rnose = 1.5 in. 
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Figure 9.116. Heat transfer calculations for 42.75° hyperboloid, R no3e = 4.489 ft. 

Case 3 

In the past, hyperboloids have been used to model the windward symmetry plane 
of the Shuttle at angle of attack (see Ref. [74], for example). As a final case, the present 
method is used to calculate the Mach 15.7 flow over a 42.75° hyperboloid, which 
approximates the windward symmetry plane of the Shuttle at an angle of attack of 
42°. The nose radius is R nose — 4.489 ft and the solution is computed for a body length 
°f Sb end — 25. The freestream conditions for an altitude of 60.56 km are p ^ = 0.00276 
psi and T, x, = 455°/?, with a variable wall temperature input. Solutions are calculated 
for noncatalytic and fully-catalytic wall conditions. Figure 9.116 compares the two 
heating rate distributions with Shuttle windward centerline measurements from STS- 
2 [31]. As expected these two calculations bracket the measured values. Note that 
the flight data is closest to the fully-catalytic result. Results presented in Ref. [74] 
also show the flight data falls closer to the fully-catalytic solution. Those results are 
for a hyperboloid geometry as well, and the solution is calculated from VSL method 
of Ref. [8]. 
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10 Conclusions 


A new approximate VSL approach to solving hypersonic flowfields about axisym- 
metl bLt bodies has been developed. The method is applied to sphe™ mrd 

hyperboloids, over 's«fto be smooth and accurate for 

perfet gH abreacting Hows. This is an improvement over the results «f Granted, 
and can be at least partially attributed to differences m the transformation van 

f01 LTMlnrpressnre relation in lieu of numerically integrating the norm^nce 
mentum equation is shown to give accurate pressure profiles outside of the pressur 
overexpansion/recompression region. As a result, predicted body pressures are g 
erallv within five percent of the VSL values. Not surprisingly, computing the norm 
component of velocity from the continuity equation, instead of using the W oxl “ a e 
l L nf Ref [291 is shown to yield profiles which are more consistent with a full 
VSL "in io^ SuL to component is small relative to the tangential velocity, 
i “f the Approximate normal velocity expression on the tangential velocity 

excellent agr.me. 

hetw^ the VSL and approximate approaches is seen d 

nress ure overexpansion/recompression region, flu* Grantz [29] and present methods 
V eld comparable results which deviate from the VSL solution. Further downstream 
L the sharp-cone limit is approached, the present technique shows good agreement 
with the VSL solution, while the results of Ref. [29] are higher. This deviation appea 
to be due to the normal velocity inaccuracies (mentioned above) having an effect on 
11 iValnv nrofiles For reacting flows, both equilibrium and nonequi- 

;';:,: n rrr,titm p L p prt,;t ^ ^ f h 

l -t L not surprising to see that the skin friction results from ™ 

approaches are in good agreement in those regmns where the heatmg agr 

Wel In the present method, a limiting form of the governing equations Wj . be ^obtained 

f - - Si? 

[luted with the present method approach their stagnation values more smoothly 
Aparins' bTtw^ He [mtnl aud VSL [9, 14) approaches for perfect 

raw - — i- •- «“ ”™ 
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the VSL algorithm. Based on the accuracy of the computed surface properties, this 
new approach could be useful in the preliminary design environment. Alternately, 
it could be used to generate an initial shock shape for more exact methods which 
require starting solutions. 
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1 1 Recommendat ions 

For this new approach to be most useful in the design environment, its extension to 
three-dimensional flows is required. The approximate approach of Ref. [29] does not 
provide the desired accuracy, so the transverse momentum equation (or an approxi- 
mate form of it) probably would need to be included in the governing equation set. 
In addition, certain facets of the ongoing research by Riley and DeJarnette [69] con 
cerning the inviscid solution of three-dimensional flows might be extended to viscous 

° W Run-time comparisons show the present technique is significantly faster than the 
VSL approach. However, especially for nonequilibrium flows, it is unclear whether this 
is due primarily to Maslen’s approximate relation or to other techniques developed 
herein to increase computational efficiency. It would be interesting to see if the run 
time would be increased by a large amount with the normal momentum equation 
included. In this scenario, the full VSL equations would be solved m a shock-normal 

^ The present method employs a marching procedure without global iterations in 
the region where the outer portion of the shock layer is supersonic. As a result, 
the physically correct phenomenon of information propagating upstream through the 
subsonic portion of the boundary layer is neglected. With the VSL approach globa 
passes followed by shock-shape smoothing allow information to propagate upstream. 
An investigation into the effect of this omission in the current approach is warranted. 

Currently, calculations cannot be made for very blunt bodies (greater than ap- 
proximately 45°). This limitation prevents the calculation of flowfields over many 
ASTV-type vehicles. The shock shape equation for the subsonic-transonic region 
would require modification before such computations could be performed. 

As a final comment, a more sophisticated procedure for determining the stream- 
wise step size in the marching region would yield significant improvements ra overall 
computational requirements. Such improvements would be most noticeable in the 
pressure overexpansion/recompression region. 


138 


REFERENCES 


References 

[1] J. D. Anderson, Jr. “A Survey of Modern Research in Hypersonic Aerodynam- 
ics”. AIAA Paper 84-1578, June 1984. 

[2] O. M. McWherter, R. W. Noack, and W. L Oberkampf. “Evaluation of In- 
viscid/Boundary Layer and Parabolized Navier-Stokes Solutions for Design of 
Reentry Vehicles”. AIAA Paper 84-0486, January 1984. 

[3] P. A. Gnoffo. An Upwind-Biased Point-Implicit Relaxation Algorithm for Vis- 
cous, Compressible Perfect-Gas Flows”. NASA TP 2953, June 1990. 

[4] W. S. Helliwell, R. P. Dickson, and S. C. Lubard. “Viscous Flow over Arbitrary 
Geometries at High Angles of Attack”. AIAA Journal, Volume 19: pages 191- 
197, February 1981. 

[5] P. A. Gnoffo. “Hypersonic Flows Over Biconics Using a Variable- Effective- 
Gamma, Parabolized Navier Stokes Code”. AIAA Paper 83-1666, July 1983. 

[6] S. L. Lawrence, D. S. Chaussee, and J. C. Tannehill. “Application of an Upwind 
Algorithm to the Three-Dimensional Parabolized Navier-Stokes Equations”. AIAA 
Paper 87-1112-CP, June 1987. 

[7] R. J . Davis. “Numerical Solution of the Hypersonic Viscous Shock Layer Equa- 
tions”. AIAA Journal , Volume 8(Number 5): pages 843-851, May 1970. 

[8] J. N. Moss. ‘Reacting Viscous Shock-Layer Solutions with Multicomponent 
Diffusion and Mass Injection”. NASA TR R-41.1, June 1974. 

[9] A. L. Murray and C. H. Lewis. “Hypersonic Three-Dimensional Viscous Shock 
Layer Flows over Blunt Bodies”. AIAA Journal , Volume 16: pages 1279-1286, 
December 1978. 

[10] R. N. Gupta, P. A. Gnoffo, and R. W. MacCormack. “Viscous Shock-Layer 
Flowfield Analysis by an Explicit-Implicit Method”. AIAA Journal , Volume 23: 
pages 723-732, May 1985. 

[11] E. C. Anderson and J. N. Moss. “Viscous-Shock-Layer Solution for Turbulent 
Flow of Radiation Gas Mixtures in Chemcal Equilibrium”. NASA TM-X 72764, 
August 1975. 

[12] J. V. Rakich, D. A. Stewart, and M. J. Lan franco. “Results of a Flight Experi- 
ment on the Catalytic Efficiency of the Space Shuttle Heat Shield”. AIAA Paper 
82-0944, June 1982. 



REFERENCES 


139 


[13] E. W. Miner and C. H. Lewis. “Hypersonic Ionizing Air Viscous Shock-Layer 
Flows Over Nonanaiytic Blunt Bodies”. NASA CR 2550, May 1975. 

[14] S. Swarninathan, M. D. Kim, and C. 11. Lewis. “Three-Dimensional Nonequi- 
librium Viscous Shock-Layer Flows Over Complex Geometries”. AIAA Paper 
83-0212, January 1983. 

[15] F. G. Blottner. “Viscous Shock Layer at the Stagnation Point with Nonequilib- 
rium Air Chemistry”. AIAA Journal , Volume 7( Number 12): pages 2281-2288, 
December 1969. 

[16] R N Gupta and A. L. Simmonds. “Hypersonic Low-Density Solutions of the 
Navier-Stokes Equations with Chemical Nonequilibrium and Multicomponent 
Surface Slip”. AIAA Paper 86-1349, .June 1986. 

[17] B. A. Bhutta and C. II. Lewis. “Recent Improvements in the Nonequilibrium 
VSL Scheme for Hypersonic Blunt- Body Flows”. AIAA Paper 91-0469, January 

1991. 

[18] S. H. Maslen. Inviscid hypersonic flow past smooth symmetric bodies. AIAA 
Journal, VolumeS: pages 1055-1061, June 1964. 

[19] W. L. Grose. An Approximate Solution to the Non equilibrium Flow in the Invis- 
cid Shock Layer about a Vehicle in Hypersonic Flight in an Arbitrary Atmosphere. 
PhD thesis, Virginia Polytechnic Institute, June 1969. 

[20] W. L. Grose. “A Thin-Shock- Layer Solution for Nonequilibrium, Inviscid Hyper- 
sonic Flows in Earth, Martian, and Venusian Atmospheres”. NASA TN D-6529, 
December 1971. 

[21] S. II. Maslen. “Asymmetric Hypersonic Flow”. NASA CR 2133, September 
1971. 

[22] E. V. Zoby and R. A. Graves, Jr. “A Computer Program for Calculating the 
Perfect Gas Inviscid Flow Field About Blunt Axisymmetric Bodies at an Angle 
of Attack of 0 Degrees”. NASA TM-X 2843, December 1973. 

[23] S K. Jackson. The Viscous-Inviscid Hypersonic Flow of a Perfect Gas over 
Smooth Symmetric Bodies. PhD thesis, University of Colorado, August 1966. 

[241 E. V. Zoby and A. L. Simmonds. “Engineering Flowfield Method with Angle of 
Attack Applications,”. Journal of Spacecraft and Rockets , Volume 22(Number 
4): pages 398-404, July-August 1985. 

1251 E V Zoby “Approximate Heating Analysis for the Windward Plane of Shuttle- 

Like Bodies at Large angles of Attack”. In T. E. Horton, editor, Thermophys.cs 


140 


REFERENCES 


of Atmospheric Entry; Progress in Astronautics and Aeronautics, volume Volume 
82, pages 229-247. AIAA, New York, 1982. 

[26] F. R. DeJarnette, II. H. Hamilton II, K. J. Weilmuenster, and F. M. Cheat- 
wood. “A Review of Some Approximate Methods Used in Aerodynamic Heating 
Analyses”. AIAA Paper 85-0906, June 1985. 

[27] A. C. Grantz. An Approximate Viscous Shock Layer Method for Calculating the 
Hypersonic Flow over Blunt-Nosed Bodies. PhD thesis, North Carolina State 
University, 1989. 

[28] A. C. Grantz, F. R. DeJarnette, and R. A. Thompson. “An Approximate Vis- 
cous Shock Layer Method for Calculating the Hypersonic Flow over Blunt-Nosed 
Bodies”. AIAA Paper 89-1695, June 1989. 

[29] A. C. Grantz, F. R. DeJarnette, and R. A. Thompson. “Approximate Viscous 
Shock-Layer Method for Hypersonic Flow over Blunt-Nosed Bodies”. Journal 
of Spacecraft and Rockets, Volume 27(Number 6): pages 597-605, November- 
December 1990. 

[30] F. M. Cheatwood and F. R. DeJarnette. “An Approximate Viscous Shock Layer 
Approach to Calculating Hypersonic Flows about Blunt-Nosed Bodies”. AIAA 
Paper 91-1348, June 1991. 

[31] D. A. Throckmorton. “Benchmark Determination of Shuttle Orbiter Entry Aero- 
dynamic Heat-Transfer Data”. Journal of Spacecraft and Rockets , Volume 20: 
pages 219-224, May- June 1983. 

[32] F. M. White. Viscous Fluid Flow. McG raw- Hill Book Company, New York, 
1974. 

[33] E. W. Miner, E. C. Andei son, and C. II. Lewis. “A Computer Program for 
Two-Dimensional and Axisymmetric Nonreacling Perfect Gas and Equilibrium 
Chemically Reacting Laminar, Transitional, and/or Turbulent Boundary Layer 
Flows”. Technical Report VPI-E-71-8, Virginia Polytechnical Institute, May 
1971. 

[34] N. B. Cohen. “Correlation Formulas and 'Tables of Density and Some Thermo- 
dynamic Properties of Equilibrium Dissociating Air for Use in Solutions of the 
Boundary-Layer Equations”. NASA TN D-194, February 1960. 

[35] C. F. Hansen. “Approximations for the Thermodynamic and Transport Proper- 
ties of High Temperature Air”. NASA TR R-50, 1959. 

[36] C. H. Lewis. “VSLNOSE: A Three-Dimensional Viscous Shock-Layer Code for 
Arbitrary Nosetips”. Technical Report Volume I: Engineering Analysis, January 
1983. 



REFERENCES 


141 


[37] W. E. Moeckel and K. C. Weston. “Composition and Thermodynamic Properties 
of Air in Chemical Equilibrium”. NACA TN 4265, 1958. 

[38] W. L. Bade. “Simple Analytical Approximation to the Equation of State of 
Dissociated Air”. ARS Journal , Volume 29( Number 4), April 1959. 

[39] S. Srinivasan, J. C. Tannehill, and K. J. Weilmuenster. “Simplified Curve Fits 
for the Thermodynamic Properties of Equilibrium Air”. NASA RP 1181, August 

1987. 

[40] S. Srinivasan and J. C. Tannehill. “Simplified Curve Fits for the Transport 
Properties of Equilibrium Air”. NASA C'R 178411, December 1987. 

[41] E. C. Anderson and J. N. Moss. “Numerical Solution of the Hypersonic Viscous- 
Shock-Layer Equations for Laminar, Transitional, and Turbulent Flows of a Per- 
fect Gas Over Blunt Axially Symmetric Bodies”. NASA TN D-7865, February 

1975. 

[42] M. C. Friedcrs, , and C. H. Lewis. “Effects of Mass Transfer into Laminar and 
Turbulent Boundary Layers over Cones at Angle of Attack”. Technical Report 
VPI-AERO-031, Virginia Polytechnical Institute, March 1975. 

[43] R. R. Thareja, K. Y. Szema, and C. H. Lewis. “Chemical Equilibrium Laminar or 
Turbulent Three-Dimensional Viscous Shock-Layer Flows”. Journal of Spacecraft 
and Rockets, Volume 20: pages 454 -460, September-October 1983. 

[44] R. N. Gupta, K. P. Lee, E. V. Zoby, J. N. Moss, and R. A. Thompson. “Hy- 
personic Viscous Shock-Layer Solutions over Long Slender Bodies-Part I. High 
Reynolds Number Flows”. Journal of Spacecraft and Rockets, Volume 27( Num- 
ber 2): pages 175 184, March-April 1990. 

[45] T. Cebeci and A. M. 0. Smith. “A Finite-Difference Solution to the Incompress- 
ible Turbulent Boundary Layer Equations by an Eddy Viscosity Concept”. In 
Computation of Turbulent Boundary Layers, Volume 1, pages 346-355. AFOSR- 
IFP-Stanford Conference, Stanford University, 1968. 

[46] T. Cebeci. “Behavior of Turbulent Flow near a Porous Wall with Pressure Gra- 
dient”. AIAA Journal, Volume 8: pages 2152-2156, December 1970. 

[47] T. Cebeci and A. M. 0. Smith. “A Finite-Difference Method for Calculating 
Compressible Laminar and Turbulent. Boundary Layers”. Journal of Basic En- 
gineering, pages 523-535, September 1970. 

[48] B. S. Baldwin and H. Lomax. “Thin Layer Approximation and Algebraic Model 
for Separated Turbulent Flows”. AIAA Paper 78-0257, January 1978. 


142 


REFERENCES 


[49] L. A. King and D. A. Johnson. “Separated Transonic Airfoil Flow Calculations 
with a Nonequilibrium Turbulence Model”. NASA TM 86830, 1985. 

[50] T. J. Coakley. “Turbulence Modeling Methods for the Compressible Navier- 
Stokes Equations”. AIAA Paper 83-1693, July 1983. 

[51] E. R. Van Driest. “On Turbulent Flow Near a Wall”. Journal of Aeronautical 
Sciences, Volume 23(Number 11): pages 1007-1011, November 1956. 

[52] F. H. Clauser. “The Turbulent Boundary Layer”. In H. L. Dryden and T. von 
Karmen, editors, Vol. IV of Advances in Applied Mechanics, pages 1-15. Aca- 
demic Press, New York, 1956. 

[53] P. S. Klebanoff. “Characteristics of Turbulence in a Boundary Layer with Zero 
Pressure Gradient”. N AC A Report 1247, 1955. 

[54] S. Dhawan and R. Narashima. “Some Properties of Boundary Layer Flow during 
Transition from Laminar to Turbulent Motion”. Journal of Fluid Mechanics, 
Volume l(Part 4): pages 418 436, January 1958. 

[55] P. S. Klebanoff. “Characteristics of Turbulence in a Boundary Layer with Zero 
Pressure Gradient”. NACA TN 3178, July 1954. 

[56] R. A. Thompson, E. V. Zoby, K. E. Wurster, and P. A. Gnoffo. “An Aerother- 
modynamic Study of Slender Conical Vehicles”. AIAA Paper 87-1475, June 
1987. 

[57] J. E. Harris. “Numerical Solutions of the Equations for Compressible Laminar, 
Transitional and Turbulent Boundary Layers and Comparisons with Experimen- 
tal Data”. NASA TR 368, 1971. 

[58] F. G. Blottner. “Finite-Difference Methods of Solution of the Boundary-Layer 
Equations”. AIAA Journal , Volume 8(Number 2): pages 193-206, February 
1970. 

[59] R. T. Davis. “Hypersonic Flow of a Chemically Reacting Binary Mixture Past 
a Blunt Body”. AIAA Paper 70-805, July 1970. 

[60] R. A. Thompson, K. P. Lee, and R. N. Gupta. “Computer Codes for the Eval- 
uation of Thermodynamic Properties, Transport Properties, and Equilibrium 
Constants of an 11-Species Air Model”. NASA TM 102602, February 1990. 

[61] R. N. Gupta, J. M. Yos, R. A. Thompson, and K. P. Lee. “A Review of Reaction 
Rates and Thermodynamic and Transport Properties for an 11-Species Air Model 
for Chemical and Thermal Nonequilibrium Calculations to 30,000 K”. NASA RP 
1232, August 1990. 



REFERENCES 


143 


[62] B. F. Armaly and K. Sutton. “Viscosity of Multicomponent Partially Ionized 
Gas Mixtures”. AIAA Paper 80-1495, July 1980. 

[631 E. A. Mason and S. C. Saxena. “Approximate Formula for the Thermal Con- 
ductivity of Gas Mixtures”. The Physics of Fluids , Volume 3(Number 5). pages 
361-369, September- October 1958. 

[641 M. D. Van Dyke and H. D. Gordon. “Supersonic Flow Past a Family of Blunt 
Axisymmetric Bodies”. NASA IRR-1, 1959. 

[651 F. G. Blottner. “Introduction to Computational Techniques for Boundary Lay- 
ers”. Sandia Laboratories Release SANl) 79-0893, September 1979. 

[66] R. D. Richtmyer. Difference Methods for Initial- Value Problems. Interscience 
Publishers, Inc., New York, 1957. 

1671 C J Riley “An Approximate Method for Calculating Inviscid Flowfields Over 

11 Three-Dimensional ^Blunt-Nosed Bodies in Hyper»nic Flow”. Master’s thes.s, 
North Carolina State University, May 1988. 

[681 C J Riley and F. R. DeJarnette. “An Approximate Method for Calculating 

11 Hypersonic Three-Dimensional Inviscid Flow Fields”. NASA TP 3018, August 

1990. 

[691 C. J. Riley and F. R. DeJarnette. “Engineering Calculations of Three-Dimen- 
sional Inviscid Hypersonic Flow Fields”. AIAA Paper 91-0701, January 1991. 

[701 R. A. Thompson and E. V. Zoby. “Flowfield and Vehicle Parameter Influence on 
Hypersonic Heat Transfer and Drag”. Journal of Spacecraft and Rockets , Volume 
27(Number 4): pages 361-368, July-August 1990. 

[71] J. W. Cleary. “Effects of Angle of Attack and Bluntness on ' 

Rate Distributions on a 15 Degree Cone at a Mach Number of 10.6 . NASA TN 

5450, October 1969. 

[79] P C Stainback, C. B. Johnson, L. B. Boney, and K. C. Wicker. “Comparison 
“ of Theoretical Predictions and Heat-Transfer Measurements for a Flight Exper- 
iment at Mach 20 (Reentry F)”. NASA TM-X 2560, 1972. 

1731 R A Thompson, E. V. Zoby, K. E. Wurster, and P. A. Gnoffo. “Aerothermody- 
namic Study of Slender Conical Vehicles". Journal of Thermophysics and Heat 
Transfer , Volume 3(Number 4): pages 361-367, October 1989. 

(741 I I Shinn J N. Moss, and A. L. Simmonds. “Viscous-Shock-Layer Heating 
Analysis for the Shuttle Windward Plane with Surface Finite Catalytic Recom- 
bination Rates”. AIAA Paper 82-0842, June 1982. 


144 


REFERENCES 


[75] D. A. Anderson, J. C. Tannehill, and R. U. Plctcher. Computational Fluid 
Mechanics and Heat Transfer. Hemisphere Publishing Corporation, New York 
1984. 

[76] J. D. Anderson, Jr. Modern Compressible Flow: with Historical Perspective. 
McGraw-Hill Book Company, New York, 1982. 

[77] J. D. Anderson, Jr. Hypersonic and High Temperature Gas Dynamics. McGraw- 
Hill Book Company, New York, 1989. 

[78] W. G. Vincenti and C. H. Kruger, Jr. Introduction to Physical Gas Dynamics. 
John Wiley and Sons, Inc., 1965. 



145 


A Reference Quantities 

The governing equations are nondimensionalized according to the following rela- 
tions: 
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and n* TtJ is the coefficient of viscosity evaluated at T* tJ . 

In addition, the following dimensionless parameters appear in the equations: 
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B Fluid Equations 


The development of the Navier-Stokes equations can be found in a variety of 
sources within the literature, so in the interest of brevity, this derivation is no . re- 
peated here. In particular, Ref. [75] presents the Navier-Stokes equations written in 
an orthogonal curvilinear coordinate system. Restricting ourselves to axisymmetric 
flow, these equations are written below in dimensional form. 
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Using the definitions of Appendix A, Eqs. (B.l) through (B.5) are rewritten below 
in nondimensional form. 
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By assuming that v and n are of the order of t and neglecting terms of higher order 

than e, the standard viscous shock-layer equations are obtained: 

continuity: 

(puhs) + -zr~ (pvhjhs) = 0 (B.ll) 
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Note: In order to keep the equations parabolic in nature, the right-hand 

side of the n-momentum equation is only of order 1. 
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C, MASLEN’S METHOD 


C Maslen’s Method 


Maslen’s pressure relation is an approximate solution to the normal momentum 
equation. Briefly, some terms are observed to be negligible in certain regions, while 
others are assumed to deviate little from their shock values. The resulting closed-form 
expression for p is a function of shock properties and stream function. The details of 
this derivation are given below, since some of the intermediate results are used in the 
present method. 

Recall from the continuity equation, Eq. (B.ll), that 

-(puh 3 ) + — (pvh 163) = 0 

This expression is satisfied by the stream function 'P, where 
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Maslen works in a(£, p) coordinate system, where 
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From Eq. (B.13), the normal momentum equation is 

( u dv dv u 2 dh\ ^ L dp^ _ n 
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Applying the transformations of Eq. (C.3) to Eq. (C.5), substituting 
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Apply Eq. (C.3) to n to get 
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By substituting this expression into Eq. (C.7) and evaluating the remaining dependent 
variables at the shock, the following relation for v is obtained 
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Differentiate with respect to p to get 
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Substitute this expression into Eq. (C.6) and evaluate the remaining dependent vari- 
ables at the shock so that 
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where the partial derivative of v with respect to £ has been neglected. Integrate this 
expression to get Maslen’s second-order pressure equation: 
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This result is essentially a truncated Taylor’s series expansion of p about its shock 
value: 

p(t,v) = MO + MOfo-i] + paCO— (c.i3) 

where 
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The streamwise derivative of p along lines of constant 7 / n appears in the governing 
equations. Differentiating Eq. (C.12) with respect to £ gives an expression for the 
streamwise derivative along lines of constant //: 
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Eq. (C.15) is the partial derivative of p along lines of constant p. However, the desired 
quantity is the partial derivative of p along lines of constant p n . Recall from Eq. (C.3) 
that 
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and from Eq. (2.1.7) 
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Equate these expressions to get 
dp _ dp dp dp 


It turns out that the expression for dp/ds, Eq. (C.4), is not well-behaved near the 
stagnation line. By making use of Eq. (C.7), this expression may be rewritten as 

dp = pub dn (C.16) 

d .< vK. df 


where dn/d£ is defined in Eq. (C.9). Equation Eq. (C.16) is well-behaved near the 
stagnation line so this form is used everywhere. 
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D Shock Layer Thickness 

With each iteration in the solution of the governing equations, a calculated value 
of the shock layer thickness, nj, may be determined from the continuity equation. As 
explained in Chapters 7 and 8, this thickness is compared with the geometric thickness 
given by the shock and body equations. When these two values are in agreement, 
the solution is advanced to the next station. This appendix provides the expression 
for determining this calculated shock layer thickness. Inherently, this expression also 
provides the relationship between 77 and rj n . 

The continuity equation, Eq. (2.1.10), can be rewritten as 


-57 ( n b p 3 pu 3 uh 3 ) - p 3 u 3 —j~ ~x~ [(t?„ - 1) puh 3 ] - p 3 (/w/Ms) = 0 (D.l) 

Integrate across the shock layer to obtain 
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Substitute for h 3 to get 
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J 0 

This yields the following quadratic for n», the normal distance from the shock to the 
body: 
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An alternate approach for calculating the shock-layer thickness is given below. 
For a given station (where s is a constant), utilizing Eq. ( . ) gives 
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which is equivalent to Eq. (D.8). Eq. (I).9) defines the relationship between rj and 
, /n . As a result, it must be used in conjunction with Maslen s pressure equation 

order to find p for a given value of ({,’/»)■ 
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E. GEOMETRIC LIMITS 


E Geometric Limits 

When the governing equations are applied on the stagnation line, their limiting 
forms (as £ ► 0) are required. These limiting expressions contain fractions in inde- 

terminant form which must be evaluated through the use of l’Hopital’s rule. This 
appendix provides the finite, nonzero values of these quantities. 

As £ — > 0 , 
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In addition to these fundamental quantities, two other terms merit attention. The 
first appears in Eq. (2.4.12), the limiting form of the streamwise pressure derivative: 
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A second term to be evaluated appears in Eq. (2.4.16), the limiting form of the 
streamwise derivative of n: 


1 dn 

H \7. Ti 


) v Sn f cos F a r,Kj 1 

^ — ) 


(E.7) 


Note that 


lim 
* 


so that 




I 1 dn 


6 ti, 0 df J 


(E.8) 
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F Shock Properties 

Boundary conditions along the body and at the shock wave (in the form of jump 
conditions) are required for the solution of the governing equations. The jump con- 
ditions across the shock vary according to the chemical nature of the flowfield. This 
appendix provides the jump conditions for perfect gas, equilibrium, and nonequilib- 
rium flows. 

In the shock-normal system 


v 


* 

CO 


VC sin I\ 


From the conservation equation 


U *oo = < = 1C COS r. 


P oo ^ OO Ps ^3 

From the momentum equation 

p* + p* u* 2 = p* -f — p* 0 k V* V* 

r co 1 roo co rs 1 Is 3 rs ' too co s 

From the energy equation 

^★2 *2 

C + -f = V + y 

Writing these equations in nondimensional form (see Appendix A) yields 


u„ = cos r. 

(F.l) 

Ps V 3 = sin r s 

(F.2) 

p s + v s sin F a = poo + sin 2 F, 

(F.3) 

e 2 sin 2 

h s 

(F.4) 

Eq. (F.3) gives 


Ps = Poo + sin 2 T s ( 1 - — ) 

V Ps) 

(F.5) 


Substituting Eq. (F.2) into Eq. (F.4) gives 


hs — h<x> T 


sin 2 r, / 1\ 

2 V ft) 


(F.6) 


The derivative of the u a with respect to ( appears in the streamwise momentum 
equation, while Maslen’s relation for the streamwise pressure gradient requires the 
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^-derivatives of p, and v a . First, differentiate Eqs. (F.l), (F.2), (F.3), and (F.4) with 
respect to £ to obtain 

(F.7) 
(F. 8 ) 
(F.9) 


du, . p 

( . — K s Sill 1 5 

dt 


dv s 


Pa 


K, COS F a + V t 


dps 


dv. 


and 


—El - = — 2 k s sin T s cos F s + v s k s cos F a sin T 3 

Ph . dv, 

— - = —k s sin F s cos I , — v s~77 
d£ “s 


(F.10) 


rpsnectivelV- • 

The closure of the above two sets of equations requires the use of the equation 

of state and the definition of enthalpy. It is advantageous to use different forms of 
these two expressions for various flowfiekl regimes. Thus, at this point the discussion 
branches into three such areas: perfect, thermally perfect, and chemically reacting 


gases. 


F.l Perfect Gas 

For a perfect gas, the specific heats are constant so that enthalpy can be expressed 


as 



/ * *rp+ 

11 s — Gp s l s 


which gives 

h s = C P ,T S = T, 

(F.ll) 

From the equation of state 

ps =: p s n at rT s 

(F.12) 

where 

7 ?*. 

7 ? = tttr 



poo 

Eqs. (F.2), (F.3), (F. 4 ), (F.ll), and (F.12) provide five equations for five un- 
knowns. They may be solved algebraically (see [76] for details) to get 

_ 2 sin 2 F, _ 7 ~ 1 (F.13) 

s_ 7 + 1 7 (7 + l)'Moo 

_ (7+ 1) sin 2 T, 

Ps ~ 2 + (7 — 1 ) sil ‘ 2 
t - E± 7 

1 S 1 

/j, 7 - I 


P. 


(F.14) 

(F.15) 
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sin I’. 


(F. 16) 


Further, since differentiating Eq. (F.14) with respect to £ gives 

dp s —4 K s p] cos r s 


dt ( 7 + l) Ml sin 3 F, 
equations Eqs. (F. 8 ) and (F.9) become 


dv * _ r 
|» ^3 COS 1 5 

d£ 


1 


dp, 

dt 


4 


[( 7 + 1) A/£,sin 2 r 3 p a J 

4/c. sin T. cos T. 


7 + 1 


-k. sin r. cos F 


7 + 1 


F.2 Nonequilibrium Flow 


(F.17) 


(F.18) 

(F.19) 


As noted in Chapter 5, for a chemically reacting mixture, the enthalpy is given by 

N s 

h = xy./t, 

i-\ 


so that 


N s 

dh = C Pj dT + J2 h >dci 

t~l 


For nonequilibrium flowfield calculations, it is assumed that the flow is chemically 
frozen through the shock wave. Thus, behind the shock 


dc{ = 0 


so that 
where 

In addition, 


dh, = C pa dT a 


c„ = C rik (T.) + (T.) 


(T,) + (T,) 

From the equation of state 

, r 

Ps = Ps-TJ-ls 

Ad n if 


7 ? = — ’±- 
,Vu n* 

CO 


(F. 20 ) 


(F .21) 


where 



F.3. EQUILIBRIUM FLOW 

Differentiating Eqs. (K.20) and (F.21) with respect to i gives 


dh,_ dT 
d{ " C ' Pi di 


dp 3 Hu 


dr ° , r dp* 

+ ‘W 
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(F.22) 

(F.23) 


d( Mu„ L 

Combining Eqs. (F.22) and (F.10), substituting the result into Eq. (F.23), and 
using that expression with Eq. (b-8) gives 

cos r. (».-/>. sin r.) 

dp, _ 'Ru dj ('ps 


T H — 

1 S I 


V 

c 


,2 

s 

Ps 


(F.24) 


Substituting Eq. (F.8) into Eq. (F.9) yields 


( }p± = _L / ( lEl _ 2k, cos r, (v, - sin l\) 
d£ v) \ di 


1 


(F .25) 


Equating Eqs. (F.24) and (F.25) gives 


dp, 

di 


K s COS T s | 

\ (^5 ps sin F 5 ) q 

1 cZ 

T 1 

— + — 

(v, - sin r.)| 

Mair _(T l 

n u \ ^5 Cps ) 



(F.26) 


Substituting Eq. (F.25) into Eq. (F.8) gives 


^i=2/c s cosr, (-- 1 

di \ps 


Ps 


r , 1 d P> 

k, cos i s H jr 

V , di 


(F.27) 


F.3 Equilibrium Flow 

For an equilibrium chemically reacting mixture of perfect gases, the definition of 
enthalpy is 

h, = h,(p s , T,) 

Typically, equilibrium air models must be employed in the form of “table look-up” 
or “curve fit” procedures. Since the energy equation is cast in terms of h, the mo e 

used here is of the form 

T, = T,{p„ h,) Z, = Z, (p s , h,) (F.28) 

From the equation of state ^ /p 2g\ 

ps — ps^s^airTs ^ 

In general, algebraic expressions for the derivatives of p. and v. with respect to { are 
difficult to obtain in the manner applied in the previous two regimes. e appr 
taken in this work is detailed in Section 4.3. 
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G Chemically Reacting Flows 

Hie governing equations as presented in Appendix B may be applied to equilib- 
rium flows, as will be discussed later in this appendix. For chemical nonequilibrium 
flows, expressing the energy equation in terms of T , rather than h is beneficial since 
the chemical finite rate equations are dependent on T. First, rewrite Eq. (B.14) as 


( u dh dh\ u dp dp 

P \T l Ys' VV dn)~~h i '¥s~ V ^n 


(G.l) 


_9q 2 f 

dn I ^ \h\ dn ' /i 3 dn 


i dh i dh 3 \ 

+ 1 — — I + p 


( du u dh\ 

\dn hi dn 


For perfect gets flow, the heat flux in the n-direction is 


q * = -k 


dT * 


dn* <ic 

so that in nondimensional form (see Appendix A), 

2 , W 

,l= - ek a^ = ,k 

where the subscript c denotes the energy flux due to conduction in the normal direc- 
tion. Since 


(G.2) 


(G.3) 


dh = C p dT 


this may be written as 


q = q c = —t‘ 


p dh 


Prdn (G ' 4) 

which is the representation used in Eq. (B.14). 

For chemically reacting flow, the effect of diffusion must be included. Thus, 


9* = 9 C * + <U 


(G.5) 


where q* d is the energy flux due to diffusion in the normal direction. This term can 
be expressed as 

«i = -/IX (G.6) 

where T>* m is the multicomponent diffusion coefficient, N a is the number of species 
present in the mixture, and c, is the mass fraction of species i. Define the Lewis 
number, which is the ratio of the mass diffusivity to the thermal diffusivity, as 




p*c v ;v* iri 


k* 


(G.7) 
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Now Eq. (G.6) may be rewritten as 


N » J’* / P 

1=1 °p/ an 


(G.8) 


Written in nondimensional form, 


q — <lc + qd 


(G.9) 


and 


^3 L* I 


(G.10) 


In the case of equilibrium flow, the conduction and diffusion quantities are com- 
bined through the following procedure. From Eqs. (G.4), (G.9), and (G.10), 


= -e 2 { k 


For a gas in local chemical equilibrium 

dci 


dT ^kLe i>m , dcA 

k ^ + S^T '^1 


dci = 


dT 


do; 

dr + — 

op 


dp 


(G.ll) 


(G.12) 


and 


dh = 


dh 

dT 


dT + 



dp 


where the subscripts P and r denote constant pressure and temperature, 
Assuming constant pressure in the normal direction gives 


(G.13) 


respectively. 




dci dci dT 
~dn = df~dn 

(G.14) 


dh _ dh dr _ dT 

dn dT dn P dn 

(G.15) 

and Eq. (G.ll) can 

be rewritten as 




' dT ^kLe iim dTdcA 
+ s C P , W) 

(G.16) 

or 

</ = -e 2 - 

f ,ar , ar\ ,. dT 

(G.17) 
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where the reaction conductivity (due to diffusion) is 

kLe i<m dc, 

k r = It-FT—hlv i 


^ c 

1=1 L 1 


P I 


dT 


(G-18) 


This definition of total conductivity (kr) can be used to define an equilibrium Prandtl 
number: 

Pr„ = ^ (G.19) 

Kt 

This quantity is also referred to as the total or reactive Prandtl number. Substituting 
into Eq. (G.16) gives 

2 fiCp dT 2 (i dh 




Pr eq dn 




Pi'eq dn 


(G.20) 


By combining Eqs. (G.20) and (G.l), it is seen that Eq. (B.14) can be used for both 
perfect gas and equilibrium flows (provided Pr fq is used). 

For nonequilibrium flows, substituting Eq. (G.ll) into Eq. (G.l) gives 


u dh dh'' 

h\ ds dn ; 


u dp dp 

h\ ds dn 


1 d 

\k dT ] 

\dn 

k dn 

- 


, ,dT ( \ dh x , 1 dh 3 \ , 

+ k~ ( - — b - — I + /* 


+ 


d_ 

dn 


^ kLe, im dc, 
2 ^— — «» 


L«=i Cp f 


dn \h\ dn h 3 dn 

+ 


dn 


( du u dhi \ 2 

\dn /ij dn J 

^ kLe^n dcj (}_d]H 1 

C PJ 1 dn \hi dn h 3 dn ) J 


(G.21) 


The variable h still appears on the left hand side of the above equation. Note that 
for a mixture of thermally perfect gases 


Ms 

h = h(c i ,T) = '£c i h i [T) 

t = l 


Therefore, 

N. N, N, .. N, 

dh = Y,c,dh, + J^hidc, = JTWdr + Y,h,dc, 

i=l i=l i— 1 i=l 

Thus, 

dh dT ^ dc, 

Since the frozen specific heat is given by 


Ns 


Cpj — 

i = 1 
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this can be written as 


dh 

ds 


dr 


N, 


= c..^ + 5> 


dc, 


P} ds ^ t'l' 1 ' ds 


Similarly, 

«A_C *L & to 

an ~ dn + t, >,n 

Substituting into Eq. (G.21) gives 


I U 

P T 


M 


' dT dc , 

C Dl — + =7 

»= 1 


' p / ds 


+ v 


_ dT ac, 

c_ + X>> 5n 

t=i 


yp / an 


u dp dp _ 

/*i as u an 




_a 

<9n 


k d Z 

dn 


+ k— 


ar / 1 a/tj i dh 


dn l /r i dn /<3 dn 


f du u dh \\ 2 

+ ^ \an /?i dn ) 


(G.22) 


_a_ 

an 


V** kLe i>m dc, 

^~c~ hi dTi 

l=\ L PJ u 


kL€i m , ac, / 1 dill 1 dh 3 

+ E-?r-^ A «— I — — + 


6' ‘an V*! an ^ * 3 an 

1=1 rj \ 

If the assumption of binary diffusion is made, 


u 


ar, A. a c , 
c ’’/ar + S ‘ar 


+ V 


,, ar aci 


)u ap ap _ 
/»! as u an 


_a_ 

dn 

+ 


*«: 

an 


+ ^ 


ar / 1 a/ii i dh 


dn \h\ dn /i 3 dn 


d_ 

dn 


N 


vA kLei2 , ac, 




A- L/C\2 j 

+ 2H-n- h ' 

i= l * 


/an n dhi 

^ \an /ti an 

aci / 1 dh\ 1 dh 3 \ 1 
‘ On \/ii dn ^ h 3 dn J J 


(G.23) 


T dn 1 ti C PJ ‘ anj ' C p/ - an \hi dn h 3 dn J J 
Based on the binary mass diffusivity, V\ 2 , the definition of the binary Lewis number 


is 


,,*c’ *T>* 
P L Pj L 12 


ljCu = (G.24) 

In addition to modifications to the energy equation, calculating flowfields with 
finite-rate chemistry requires the solution of the species conservation equation for 
each chemical species in the mixture. An excellent derivation of this equation is 
presented in [77]. The final result is repeated here: 


f u* dci * ac, 

I U *77 

.h x ds* dn* t 


- tb * + (G.25) 

" i<; * On* J ' \hidn*^ h* 3 dn*) V 

where w* is the rate of production term for species i and J* is the diffusion mass flux 
term of species i. This binary diffusion mass flux term is defined by Fick’s law to be 

7* — —o*'V* —— 

Ji ~ 9 Ui2 dn* 


(G.26) 
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which can lx; written as 


j: = 


c 


, Oc-i 

V 


n" Oc t 

ErOn * 


Using the definitions of Appendix A, Eq. (CJ.25) is rewritten below in 
sional form: 


it dci ()c t 

h [ Os + On 


w, 


.JM 


On 


- i j. 


I i)h. 


+ 


1 0h-\ 


li i O n h \i O n 


whore 


, T _ k , Or, 

t'l ,v 

(•'P j On 


Substituting Eq. ((1.29) into Eq. (G.23) and rearranging yields 


(>('> 


it or i)T\ u op op 


N, 


P f 1 h\ Os 


+ v- 


On 


h \ 0* 


, + >> 


On 


{ 0 

i 

\0n 

[oTi\ 


, OT ( 1 dhx 1 Oh‘C 

+ k— I + 


On \h j On // ;i On 


t=\ 


+ /< 


a 


0c t 


Oc; 


h\ Os ^ V On 


( On u dh\ 

On li\ On 


Ms 


i=l 


] On 


+ J< 


1 Oh | 1 />//(' 

+ 


tv. 


h i 5/; /r t <9n 


C 2 YJC — 

tx Pi ^ 


Combining Eq. (G.28) and Eq. (G.30) gives 


f a OT 0T\ u Op Op 

,< n\u l TC v ^)-h l lir v Yn = 


0_ 

On 


OT 

On 


+ ic- 


rn' ( 1 Oh i 1 Oh-i 


+ 


( Ou u Oh i 


On \ li | On h ’i On 

OT 

On 


+ I 1 IT ~ 


\0n h\ On 


/V, uyr N, 


i = 1 


t=l 


(G.27) 

nondimen- 

(G.28) 

(G.29) 

(G.30) 


(G.31) 
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H Reaction Rates 

i 

In nonequilibrium flows, the chemical reactions proceed at a finite rate. The 
production terms (w t for each species i in the mixture) appear in the energy equation 
(when formulated in terms of temperature) and in the species continuity equations. 
A discussion of these production terms is piesented below. 

Consider as an example, the dissociation of diatomic oxygen: 

0 2 + M r — L 20+ M r (H.l) 

where M r is a collision partner (or catalytic third body). Catalytic bodies are those 
collision partners which do not undergo a chemical change during the reaction. From 
empirical results, the rate of formation of 0 can be written as 

10] = 2k) [0,1 \U,\ (H.2) 

where k) is the forward reaction rate coefficient which is a function of T only, and 
each bracketed term represents the concentration of that substance (in moles/cm ). 
The reverse (backward) reaction is 


02 + M r h- 20 + Mr 

(H.3) 

so that the rate of reduction of 0 is 


± [0] = -2k; [Of [M.1 

(H.4) 

Combining these two elements gives 


0 2 + M r -20+ M r 

(H.5) 

so that the net rate of 0 formation is 


4- [0] = 2 k) [0 2 ] [M r ] - 2kt [Of [M r ] 
dt * 

(H.6) 

For the more general case of a multicomponent gas with N s distinct reacting 
chemical species and N r simultaneous chemical reactions, the stoichiometric relations 
for the overall change from reactants to products are: 

N, Nt 

^ £ H,,x- 

1=1 

(H.7) 
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where r — 1,2, ..., N r and N t is equal to the total number of the reacting species ( N , ) 
and catalytic third bodies ( N c ). The quantities a^ r and /?, )r axe the stoichiometric 
coefficients for reactants and products, respectively. The variable X* denotes the 
concentration of species i (or catalytic third body i — N s ): 



=P ‘M, ’ 

moles 


(H.8) 


cm 3 


In the above example, 





ftOj = 1 

ato - 0 


<*M r ~ 1 


P<h = 0 

II 

O 


(I M r = 1 



For the r-th reaction of Eq. (II. 7), the lorward rate of production is 


' dX*\ Nc 

r) =(A>-«,>)^n(7./r r 


dt* 


while the backward rate is 


(H-9) 


'dX*\ Nc 

jr i = ~ (A? - °v) k l II (7 ip*f'- r 

i=i 


dt* 


(H.10) 


so that the net rate of production [78] is 






N c 


N c 


k ir n (7i7>*)“ ,,r - k n ( 7 ./) 




t=i 


1=1 


(H.ll) 


where p* is in gm/cm 3 . The net mass rate of production of the i-th species may be 
obtained by summing (Eq. H.ll) over all the reactions N r : 


or 



(H. 12) 


(H.13) 


«■), = n < 7 i />*)“•■' = n ( 7 ,)”-’ 

P i=l t=l 

^ £ iW)*' = </>'* ft (*)*' 

P 1=1 i= 1 


N c 

Ot T ^ ^ t r 1 
* = 1 


where 


169 


Table H.l. Chemical Reactions 


V 

Reaction 

M T 

1 

02 + Mi ?=* 20 + Mi 

N 2 , 0 2 , 0, NO, N 

2 

N 2 + M 2 ^ 2N + M 2 

N 2 , 0 2 , 0, NO 

3 

A 2 + M 3 ^ 2rV T 4/3 

N 

4 

jVO + M a ?= iv + 0 + m, 

N 2 , 02, 0, NO , AT 

5 

A"0 + 0 02 + -V 


6 

At 2 + 0 ^ iVO + JV 


7 

AT + 0 ^ NCR + e~ 

1 


N c 

1=] 

The mole-mass ratio (or mole number), 7 m 1S defined as 

_c^ = xr 

Mi P * 

N 


7.= 


, for i = 1,2, ..., N s 

iV 3 

, fori = iV. + l,-,W, 


(H.14) 


n— 1 


where Zr-w... is the catalytic efficiency of third body M- At, for specie n. Written in 
modified Arrhenius form, the reaction rates are 


fej f = A,, r T* s ''exp(-7’o / ,/T*) , 

k‘ hr = A t , r T‘ B ‘' eM-njT') . 


1 

^ cm 3 ^ 

Or" 

s 

1 K mole 


1 j 

1 cm 3 ^ 

A-' 

s ' 

^mo/e y 



(H.15) 

(H.16) 


"'TnthisstLdyfi seven-species model [15] for air is used. Table H.l lists the reactions 
and third bodies which are considered. The efficiencies of the catalytic third bodies 
are given in Table H.2 for the various reactions. Finally, the forward and backward 
rate coefficients are presented in Table H.3. 
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Table H.2. Third Body Efficiencies Relative to Argon 


Catalytic 

Bodies 

i-N, 

Efficiencies 

N- 2 

l = 1 

o 2 

i = 2 

N 

i = 3 

0 

i = 4 

wo 

i = 5 

wo+ 

i = 6 

A/, 

1 

2 

9 

1 

25 

1 

0 

m 2 

2 

2.5 

1 

0 

1 

1 

0 

M 4 

3 

1 

1 

20 

20 

20 

0 

e~ 

4 

0 

0 

0 

0 

0 

1 


Table H.3. Chemical Rate Coefficients 


r 

K Jr 

[cm 3 / mole-sec] 

Hr 

[art* /mole-sec] or [cm 6 /mole 2 -sec] 

i 

3.61 x 10 18 7 ,_1 -°exp( -5.94x10V'/ 1 ) 

*5 

O 

1 

1C 

O 

X 

o 

CO 

2 

1.92xl0 17 T-°- 5 exp(-1.131 x 10 5 /T) 

1.09x 10 l6 T _o - 5 

3 

4.15xl0 22 T ,-1 - 5 exp(-1.131 x 10 5 /T) 

2.32xl0 21 T -1-5 

4 

3.97 x 10 2O T -1 ' 5 exp( -7.56xl0 4 /T) 

1.01xl0 2o T -1 - 5 

5 

3.18xl0 9 T 10 exp( -1.97xl0 4 /T) 

9.63xl0 11 T’°- 5 exp(— 3.6 x 10 3 /T) 

6 

6.75xl0 13 exp( -3.75X10V7 1 ) 

1.50xl0 13 

7 


].80xl0 ,9 T’ _1 ° 




171 


I Species Thermodynamic 
Properties 

For nonequilibrium flows, species thermodynamic properties are required for each 
species present in the fluid. Ref. [60] has provided curve fits for C* and h* for the 
temperature range of 300 A < T* < 30000 A: 
specific heat: 


C 


Pi 


n u * 


A + A 2 T* + A 3 T *' 2 + A 4 T + AsT 


^+4 


static enthalpy: 


A 2 T * . A3T* 2 , A 4 T * 3 , A 5 T* 4 l Ae 

J 1 * 


^298.1,5 = A . + + 1^— + -V- + 

K u *T* A + 2 3 4 5 


(1.1) 


(1.2) 


The subscript . denotes the i-th species as given in Table 5.1. The temperature range 
cited has been divided into five sub-ranges with each sub-range curve fit with the 
above polynomials. Coefficients Ax through As are presented in Table 1.1, grouped 
according to these sub-ranges. The Aj coefficients are related to the heats of forma- 
tion of their respective species (and therefore invariant with temperature). They are 
located in Table 1.2, along with the species heats of formation. 

Note that the enthalpy from Eq. (1.2) is referenced to a temperature of 298. 15 A. 
The specific enthalpy for a reference temperature of 0A is desired. This can e 
calculated using the relation 





+ Ah* - A/C + 298.15 C* 


(1.3) 


where the first term on the right hand side is obtained from the above curve fits. 
The second and third terms (which are the heats of formation at 0A and 298.15 A, 
respectively) are provided in Table 1.2, while the perfect gas values of C* are presented 

in Table 5.1. 
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Table 1 . 1 . Constants for Polynomial Curve-fits of Thermodynamic Properties 


Species 

Ax 

Ao 

A 3 

Aa 

As 




300 A' < T* < 1000 A' 


N 2 

.36748x10* 

-.12081 xlO" 2 

,23240xl0- 5 

-.63218x10-° 

-.22577xl0"* 2 

Oi 

.36146x10* 

-.18598xl0 -2 

.70814x 10 -5 

-.68070x10"® 

.21628x10-** 

N 

.25031x10* 

-.21800x10-“ 

.54205x 10 -7 

-.56476x10"*° 

.20999X10 - * 3 

0 

.28236x10* 

-.89478xl0 -3 

.83060x10-° 

-.16837x10-° 

- 73205xl0"* 3 

NO 

.35887x10* 

-.12479xl0" 2 

.39786xl0" 5 

-.28651x10-® 

.63015xl0" 12 

N0 + 

.35294x10* 

-.30342 xl0“ 3 

.38544x10-° 

.10519x10-® 

-.72777X10"* 2 

e” 

.25000x10* 

.00000x10° 

.00000x10° 

.00000x10° 

.00000x10° 




1000 A' < r r < 6000 A' 


N 2 

.32125x10* 

.10137 x 10" 2 

-.30467x10-° 

.41091x10-*° 

-.20170x10-*“ 

0? 

.35949x10* 

.75213xl0 -3 

-.18732x10-° 

.27913x10"*° 

-.15774x10-*“ 

N 

.24820x10* 

.69258xl0 -4 

— .63065xl0 -7 

.18387x10" 10 

-.11747x10-*“ 

0 

.25421x10* 

-.27551x10"“ 

-.31028x10-® 

.45511x10-** 

-.43681x10-*® 

NO 

.32047x10* 

.12705x 10 -2 

-.46603 xlO" 6 

.75007x10"*° 

-.42314x10“*“ 

NO* 

.32152x10* 

.99742xl0" 3 

-.29030x10"° 

.36925x10-*° 

-.15994x10"*“ 

e~ 

.25000x10* 

.00000x10° 

.00000x10° 

.00000x10° 

.00000x10° 



6000 K < T* < 15000 A' 


No 

.31811x10* 

.89745xl0 -3 

-.20216x10-° 

.18266x10"*° 

-.50334x10-*® 

Oo 

.38599x10* 

,32510xl0" 3 

-.92131x10-® 

-.78684 xlO" 12 

.29426x10-*° 

N 

.27480x10* 

-,39090xl0" 3 

.13380x10-° 

-.11910x10-*° 

.33690x10-*® 

0 

.25460x10* 

-.59520x10"“ 

.27010xl0 -7 

-.27980x10-** 

.93800x10-*° 

NO 

.38543x10* 

.23409xl0" 3 

— .21354xl0 -7 

.16689x10-** 

-.49070x10-*° 

NO* 

.26896x10* 

.13796xl0" 2 

-.33985x10"° 

.33776x10" 10 

-.10427x10-*“ 

e~ 

.25000x10* 

.00000x10° 

.00000x10° 

.00000x10° 

.00000x10° 



15000 A' < 7** < 25000A' 


N 2 

.96377x10* 

— .25728xl0" 2 

.33020x10“° 

-.14315x10"*° 

.20333x10-*® 

o 2 

.34867x10* 

.52384 xlO' 3 

— ,39123x 10 -7 

.10094x10-** 

-.88718xl0-* 7 

N 

-.12280x10* 

.t9268xl0" 2 

-.24370x10"° 

.12193x10"*° 

-.19918x10-*® 

0 

-.97871 xlO- 2 

.12450x 10" 2 

-.16154x10-° 

.80380x10-** 

-.12624x10-*® 

NO 

.43309x10* 

-.58086x10-“ 

,28059xl0“ 7 

-.15694x10"** 

.24104x10-*° 

N0+ 

.59346x10* 

-.13178xl0" 2 

.23297x10-° 

-.11733x10"*° 

.18402x10-*® 

e~ 

.25000x10* 

.00000x10° 

.00000x10° 

.00000x10° 

.00000x10° 



25000 K < T* < 30000 A' 


Ni 

-.51681x10* 

.23337 x 10 -2 

-.12953x10"° 

.27872x10-** 

-.21360x10"*° 

o 2 

.39620x10* 

.39446x10“ 3 

— .29506x 10“ 7 

.73975x10" 12 

-.64209X10 - * 7 

N 

.15520xl0 2 

-.38858x10- 2 

.32288x10-° 

-. 96053 x 10" ** 

.95472x10" 16 

0 

. 16428 xlO 2 

-.39313xl0 -2 

.29840x10-° 

-.81613x10-** 

.75004x10-*° 

NO 

.23507x10* 

. 58643x10-° 

-,31316x 10 -7 

.60495 xlO" 12 

-.40557xl0-* 7 

NO* 

-.51595x10* 

.26290 xlO" 2 

-.16254x10-° 

.39381x10"** 

-.34311x10-*° 

€~ 

.25000x10* 

.00000x10° 

.00000x10° 

.00000x10° 

.00000x10° 




Table 1.2. Heats of Formation 


Species 

\kcal/gm-mole] 


N-2 ”1 

0.0 

0.0 

— .10430xl0 4 

0 2 

0.0 

0.0 

-. 10440 xlO 4 

N 

112.973 

112.529 

,56130 x10 s 

0 

59.553 

58.984 

,29150 x10 s 

NO 

21.580 

21.457 

.97640 xlO 4 

NO* 

236.660 

235.180 

,11840xl0 6 

e~ 

0.0 

0.0 

-. 74542 xlO 3 



174 


J. SPECIES TRANSPORT PROPERTIES 


J Species Transport Properties 


For nonequilibrium flows, species transport properties are required for each species 
present in the fluid. For the temperature range of 1000A < T* < 30000R", the species 
transport properties are defined by the curve fits of Ref. [60]. The expression for 
viscosity is 



/ji* 


/L. In T 


+ < 


[gm/cm-s 


(j.i) 


where the species coefficients are given in Table J.I. The frozen thermal conductivity 
is given by 



r [A ki (In T*) 3 + B k . (In T*) 2 + C\ In T* + D k \ 


[ cal/cm-s-K ] (J.2) 


with the species coefficients located in Table J.2. 

The ionic species viscosities and frozen thermal conductivities as defined by the 
tabulated coefficients are for the limiting electron pressure, which is defined as 

P*t m a: = - 0975 (^jjo) > [«H ( J - 3 ) 

These values should be corrected by employing the following formula [61]: 

Pi (pW _ h (p;,) 

P> {P*lmar) k ' (Pc/m.r 
The local electron pressure is 

P*i = N € K b T* , [atm] (J.5) 



Table J.I. Constants for Viscosity Curve-Fits 


Species 

A, 


c. 

N-2 

0.0203 

0.4329 

-11.8153 

02 

0.0484 

-0.1455 

-8.9231 

N 

0.0120 

0.5930 

-12.3805 

0 

0.0205 

0.4257 

-11.5803 

NO 

0.0452 

-0.0609 

-9.4596 

NO + 

0.0 

2.5 

-32.0453 

e~ 

0.0 

2.5 

-37.4475 
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Table J.2. Constants for Frozen Thermal Conductivity Curve-Fits 


Species 

Ak 

B k 

Ck 

D k 

E k 

iV 2 

0.03607 

-1.07503 

11.95029 

-57.90063 

93.21782 

o 2 

0.07987 

-2.58428 

31.25959 

-166.76267 

321.69820 

N 

0.0 

0.0 

0.01619 

0.55022 

-12.92190 

0 

0.0 

0.0 

0.03310 

0.22834 

-11.58116 

NO 

0.02792 

-0.87133 

10.17967 

-52.03466 

88.67060 

N0+ 

-0.06836 

2.57829 

-35.72737 

219.09215 

-519.00261 

€~ 

0.0 

0.0 

0.00032 

2.49375 

-27.89805 


where N c is the electron number density and A'b is Boltzmann’s constant In the 
present work, an electron pressure of l aim is assumed everywhere so that the above 

correction reduces to 


Hi 


m _ h 

(Pt lm ar) k ' 


2 



/ T* \ 4 rn l 

^ \ 8 / 3 l 

In 

rKiooo) +i - 52 ( 

aoooJ J 


(J.6) 
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K. CONIC SHOCK SHAPE 


K Conic Shock Shape 

Define the shock shape with a conic equation: 

r 2 — 9 n _ /, ,2 

Differentiating Eq. (K.l) with respect to £ gives 

<fr. , , cos 1\ 

-77 = (6 a - = sin I' 4 

Differentiate Eq. (K.2) with respect to £ and rearrange to get 

• 2 n 2 p 

sin 1 3 -f b s cos 1 3 


cos r 


(K.l) 


(K. 2 ) 


(K. 3 ) 


Now differentiate Eq. (K. 3 ) with respect to £ to get 

d / K s r, \ _ 2/c a sin F s cos T s (b s - 1) 
Vcos 1 

Utilizing the product rule of differentiation, 

r. \ 

+ 


-1 

( \ 

\ d 1 

f * ) 


\ cos r 3 ) 

1 - K ’d{ 1 

[cos r 3 J 


r * 


cos r* 


(K. 4 ) 


(K. 5 ) 


so that 


da, cos r, 

d£ r 4 


2« s sin I', cos r, (b 3 - 1) - k 3 (tan T, - k s sin T 4 ) j 


That is, 


< 1 k s 

~di 


= K. 


. cos 2 r s sin r 4 

2 sin T s (b 3 - 1) 


(*-- 


cos 2 r. 


Substituting Eq. (K. 3 ) and rearranging gives 


d* a o • r cos 2 r, 
— sin 1 . 


(b, ~ 1) 


d£ r 3 

On the stagnation line, Eq. (K. 7 ) has the following limiting form: 

1 da. 


(K.6) 


(K. 7 ) 


I'™ I cosT. d( 


3k,,' (6.-1) 


(K.8) 
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L Cubic Shock Shape 

Define the shock shape in the region between the previous and current marching 
stations with the following cubic equation: 


= fo + r, (Az,) + (A z,) 2 + (A: 


(L.l) 


where 


A 


s k - 1 


The axial position of the current station is - z Sh , 
Zs = z Differentiating Eq. (L.l) with respect to 


while the previous station is at 
gives the shock slope: 


= r, + fj (Az.) + ~ (Az,) 2 = tan r. 


(L.2) 


Differentiating Eq. (L.2) with respect to c 3 gives 

44 = r 2 + h (A z,) = -k s cos -3 T s 


(L.3) 


One more differentiation yields 



(L.4) 


Now focus on the previous station (fc- 1), where this current shock segment meets 
the previous segment (between stations k - 2 and k - 1). The shock geometry here 
was determined from this previous segment. Requiring a continuous position, slope, 
and curvature across this juncture (where z a = Zs^) gives 


r 0 = r, k _ { 



fc - 1 


d 2 r s 


r 2 = 


dz] 


fc - 1 


However, from Eq. (L.4), note that the curvature derivative is not continuous between 
shock segments. In fact, the derivative’s value is invariant along the length of a given 

segment . 

Eq. (L.3) can be rewritten as 


1 

d 2 r. 

d 2 r. 


A; s 

dz] 

dz 2 

fc- 1. 


(L.5) 


so that Eq. (L.4) becomes 


d 3 r a _ 1 

d 2 r a 

d 2 r s 


dz ^ Az s 

dz 2 

dz] 

fc- 1- 


(L-6) 
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L CUBIC SHOCK SHAPE 


which is simply a two-point backward difference representation of < Pr,/dz 2 . Substitute 
this definition of r 3 into Eq. (L.l) to get 


r . = r, 


Jt— i 


* dr, 

+ IT. 


+ 


Az 2 


jt-i 


2^. 


dz 2 


+ 


Jb— 1 


<Pr, 

~dz J 

3 j 


(L.7) 


which defines the shock segment in terms of values at the previous station. Similarly, 
Eq. (L.2) can be written as 


di' s 

dr. 

Az, 

d 2 r. 

d 2 r. 

dz. 

dz. 

+ 9 

k—\ 

dz] 

+ dz 2 

k - 1 UZ s J 


to yield an expression for the shock slope between stations k — 1 and k. 


(L.8) 
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